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Concurrent  MR-NIR  Imaging  for 
Breast  Cancer  Diagnosis 


4 


Birsen  Yazici 


I  INTRODUCTION 

Near  infrared  (NIR)  diffuse  optical  imaging  provides  quantitative  functional  information  from  breast  tissue  that 
can  not  be  obtained  by  conventional  radiological  methods.  NIR  techniques  can  provide  in  vivo  measurements 
of  oxygenation  and  vascularization  state,  the  uptake  and  release  of  molecular  contrast  agents  and  chromophore 
concentrations  with  high  sensitivity.  There  is  considerable  evidence  that  tumor  growth  is  dependent  on  angio¬ 
genesis  [1]-  [3],  and  that  tumor  aggressiveness  can  be  assessed  from  its  increased  number  of  new  vessels  and 
reduced  oxygenation  state  relative  to  normal  breast  tissue  and  benign  breast  lesions  [4]-  [6].  NIR  diffuse  optical 
tomographic  (DOT)  methods  has  the  potential  to  characterize  angiogenesis  related  vessel  density  as  it  measures 
the  total  hemoglobin  concentration  and  provide  the  ability  to  differentiate  between  benign  and  malignant  lesions 
based  on  oxygen  saturation.  Furthermore,  NIR  methods  are  non-ionizing,  relatively  inexpensive  and  can  be  made 
portable. 

The  diagnosis  and  management  of  cancer  involves  several  stages  where  magnetic  resonance  (MR)  plays  a 
valuable  and  growing  role.  MRI  of  the  breast  is  now  a  routine  part  of  the  clinical  care  in  many  centers  [9]- 
[11].  Magnetic  Resonance  imaging  (MRI)  is  indicated  in  patients  with  inconclusive  clinical  and/or  mammographic 
examinations.  Patients  that  may  benefit  include  women  with  radiographically  dense  breasts,  and  high  risk  potential 
population  [12]-  [13].  MRI  possesses  less  than  10%  contrast  for  soft  tissue  pathology  [14].  Gadolinium  (Gd) 
enhanced  MRI  offers  much  better  contrast  and  is  specific  for  tumor  vessel  imaging.  Flowever,  the  signal  in  the 
Gd-MRI  arises  from  the  larger  vessels  as  the  contrast  agent  is  flushed  out  of  the  vascular  bed  of  the  tumor  [15]. 
In  comparison,  NIR  measurements  of  absoiption  have  extremely  high  contrast.  It  was  reported  that  5%  change 
in  vascular  density  as  measured  histologically  in  ductal  carcinomas  leads  to  approximately  300%  contrast  in  NIR 
absoiption  coefficients  [7].  Furthermore,  there  are  studies  suggesting  that  the  kinetics  of  contrast  enhanced  optical 
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spectroscopy  provides  information  about  the  cellular  spaces  [8].  On  the  other  hand,  NIR  DOT  suffers  from  poor 
spatial  resolution  and  as  such,  it  is  unlikely  that  NIR  imaging  will  be  a  stand-alone  screening  method  in  the  general 
population.  Therefore,  we  believe  that  the  concurrent  MR  and  NIR  imaging  brings  together  the  most  advantageous 
aspects  of  the  two  imaging  modalities  (structural  and  functional).  In  the  future,  we  envision  that  this  multimodality 
imaging  approach  will  lead  to  high  resolution  hemoglobin  tomography  and  comprehensive  quantitative  functional 
tissue  characterization  to  differentiate  malignant  and  benign  tumors. 

In  this  project,  the  clinical  studies  are  performed  using  the  novel  MR-NIR  hybrid  time-resolved  spectroscopy 
(TRS)  imager  and  fast  Indocyanine  Green  (ICG)  enhanced  spectroscopic  imager  developed  by  Dr.  Chance,  a  Co-PI 
of  this  proposal,  at  the  University  of  Pennsylvania  (UPenn),  Biophysics  Department,  Diffuse  Optical  Imaging  and 
Spectroscopy  Laboratory. 

The  central  hypothesis  of  this  project  is  that  the  concurrent  MR-NIR  diffuse  optical  tomographic  methods  coupled 
with  fast  contrast  enhanced  NIR  spectroscopic  methods  provide  fundamentally  new  quantitative  functional  and 
structural  information  for  breast  cancer  tumor  characterization  and  detection.  This  new  information  can  be  obtained 
by  novel  modeling,  analysis  and  data  fusion  methods  from  the  tomographic,  temporal  and  cellular-based  contrast 
measurements,  which  exploit  fast  imaging  techniques  together  with  TRS  tomographic  methods.  In  this  project, 
we  investigate  new  methods  for  multi-modality  high  spatial  resolution  hemoglobin  tomography,  pharmacokinetic 
modeling  of  molecular  contrast  agents  based  on  fast  NIR  spectroscopy  and  analysis  of  structural  and  functional 
information  provided  by  MR  and  NIR  imaging  methods  for  breast  cancer  detection  based  on  receiver  operating 
characteristics  methodology.  Specific  aims  of  the  project  are  as  follows: 

•  Aim  1:  Utilize  a  priori  anatomical  information  provided  by  MRI,  to  reconstruct  3D  high  resolution  hemoglobin, 
water  and  lipid  concentration,  and  oxygen  saturation  images  directly  from  6  wavelength  time  resolved  optical 
measurements.  Evaluate  improvements  in  image  reconstruction  between  that  of  stand-alone  NIR  and  concurrent 
MR-NIR  measurements  using  water  and  lipid  images  obtained  from  MRI. 

•  Aim  2:  Develop  a  compartmentalized  pharmacokinetic  modeling  of  ICG,  optical  contrast  agent,  and  extract 
quantitative  parameters  that  can  characterize  tumor  metabolism  and  angiogenesis.  Compare  ICG  kinetics  with 
the  Gadolinium,  MR  contrast  agent,  kinetics  and  biopsy  findings. 

•  Aim  3:  Evaluate  accuracy  of  breast  cancer  diagnosis  based  on  the  quantitative  functional  information  extracted 
from  stand-alone  NIR  system.  This  information  includes  hemoglobin,  water  and  lipid  concentration,  optical 
scatter  power  and  oxygen  saturation  images,  and  ICG  pharmacokinetic  parameters.  Evaluate  the  added  value 
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of  ICG  kinetic  parameters  in  breast  cancer  diagnosis. 

•  Aim  4:  Combine  NIR  based  breast  cancer  diagnosis  features  with  the  systematic  MR  breast  architecture  and 
kinetics  interpretation  model  developed  by  Dr.  Nunes,  M.D,  Co-PI  of  this  proposal,  to  evaluate  the  sensitivity 
and  specificity  of  concurrent  MR-NIR  imaging  method.  Compare  results  with  that  of  stand-alone  MR  and  NIR 
results. 

In  the  following  sections,  we  will  provide  detailed  description  of  our  current  research  in  line  with  the  statement 
of  work  (SOW)  and  the  aims  outlined  above.  For  the  period  of  June  1st,  2007  to  May  31st  2008,  SOW  includes 
only  the  first  two  aims  of  the  project. 


II  BODY 

A.  AIM  2  -  Direct  Reconstruction  of  Optical  Fluorophores 

1 )  Spatially  Resolved  Compartmental  Modeling  and  the  Two-compartment  Model  for  ICG:  In  general,  a  com- 
partmental  model  is  given  by  time-dependent  coupled  ODEs  and  a  measurement  model  [1 8]— [20].  Such  a  model 
can  be  extended  to  include  spatial  variations  in  a  straightforward  manner. 

Let  O  C  M3  denote  the  domain  of  interest.  For  an  n-compartment  model,  let  C(r,t )  E  M"  represent  the 
concentration  vector  in  different  compartments  at  location  reO,  and  at  time  t  E  [To,Ti];  and  let  an(r )  denote 
the  parameter  vector  whose  elements  are  the  pharmacokinetic  rates  and  volume  fractions  at  location  r.  Then,  a 
spatially-resolved  compartmental  model  can  be  expressed  as  the  following  state-space  model: 

C(r,t)  =  K(an(r))C(r,t),  (1) 

a(r,t )  =  V{an(r))C(r,t)  te[T0,T j],  rell,  (2) 

where  C  denotes  the  elementwise  time-derivative  of  C,  JC(r,  an)  E  MnXTl  is  the  system  matrix  whose  entries  are  the 
pharmacokinetic  rates,  a(r,t )  is  the  total  fluorophore  concentration  at  location  r  and  time  t,  and  V(r,an)  E  Mn 
is  the  vector  comprised  of  volume  fractions. 

Although  our  study  is  applicable  to  pharmacokinetic  modeling  of  any  optical  fluorophore,  in  the  following 
subsection,  we  will  specifically  discuss  the  spatially-resolved  pharmacokinetic  modeling  of  ICG  using  a  two- 
compartment  model  due  to  its  relevance  to  breast  cancer  studies. 

2)  Two-compartment  Model  for  the  ICG  Pharmacokinetics:  ICG  is  an  optical  dye  commonly  used  in  retinopathy 
and  hepatic  diagnostics.  Given  its  low  toxicity  and  FDA  approval,  it  has  recently  been  utilized  as  a  blood  pooling 
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agent  for  the  detection  and  diagnosis  of  cancerous  tumors  in  conjunction  with  NIR  optical  methods  [34],  [35], 
[37],  [41].  In  normal  tissue,  ICG  acts  as  a  blood  flow  indicator  in  tight  capillaries  of  normal  vessels.  However  in 
tumors,  ICG  may  act  as  a  diffusible  (extravascular)  flow  in  the  leaky  capillary  of  cancer  vessels.  Therefore,  the 
pharmacokinetics  of  ICG  can  be  used  for  tumor  detection,  diagnosis,  and  staging.  Fig.  1  shows  the  two-compartment 
model  for  the  ICG  kinetics.  The  two  compartments  are  comprised  of  plasma  and  extracellular-extravascular  space 
(EES).  Spatially  resolved  ICG  transition  between  plasma  and  the  EES  can  be  modeled  using  the  following  coupled 
ODEs: 


1 

kout  (  ^  )  kin  (v  ) 

i - 

■40 

i _ 

3" 

e-t- 

1 _ 

kout(r)  ( kin  +  A^e/m)(r) 

CP{r ,  t) 

te[T0,T i],  rell, 


(3) 


where  kin{r )  and  kout{r)  are  the  spatially-resolved  pharmacokinetic-rates  that  govern  the  leakage  into  and  the 
drainage  out  of  the  EES,  keim(r)  describes  the  ICG  elimination  from  the  body  through  kidneys  and  liver.  The 
vector  C(r,t )  in  (1)  and  (2)  is  composed  of  Cp(r,t),  and  Ce(r,t),  representing  the  ICG  concentration  in  plasma 
and  the  EES  at  r  E  H  and  t  E  [To,  Tj],  respectively. 


Capillary 

Cp?  Vp 

kelm 

kin 

kout 

. 

EES 

Ce  ? 

Fig.  1.  Block  diagram  of  the  two-compartment  model  for  the  ICG  pharmacokinetics. 


The  ICG  concentration  in  tissue,  a(r,  t),  is  given  as  a  linear  combination  of  the  ICG  concentration  in  plasma 
and  the  EES 

a(r,t)  =  vp(r)Cp(r,t)  +  ve(r)Ce(r,t)  te[T0,Ti],  r  E  (4) 

where  vp(r),  ve(r)  are  spatially-resolved  plasma  and  EES  volume  fractions,  respectively.  Here,  the  unknowns  are 
concentrations  in  different  compartments,  pharmacokinetic  rates,  and  volume  fractions. 
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Fig.  2.  Total  fluorophore  concentration  versus  slow-time  variable,  t. 

We  combine  all  the  pharmacokinetic  parameters  into  a  single  vector  and  define 

02(f)  =  [kin(r)  kout(r)  keim(r)  vp{r)  ve(r)]T .  (5) 

3)  Concentration-to-Measurement  Map  for  FDOT:  In  this  work,  the  quantity  we  wish  to  reconstruct  is  the 
spatially  resolved  pharmacokinetic -rate  images  from  a  sequence  of  boundary  measurements  obtained  by  diffuse 
optical  tomographic  methods.  To  do  this,  we  first  develop  a  map  which  relates  sequence  of  boundary  measurements 
to  the  spatially  resolved  fluorophore  concentrations.  We  call  this  map  concentration-to-measurement  (CTM)  map. 

A  suitable  CTM  map  can  be  developed  based  on  a  photon  propagation  model  in  fluorescing  medium.  We 
use  diffusion  approximation  of  radiative  transfer  equation  to  model  photon  propagation  where  the  propagation  of 
excitation  and  emission  light  are  modeled  by  two  coupled  diffusion  equations. 

To  relate  time-varying  fluorophore  concentrations  to  a  sequence  of  boundary  measurements,  we  parameterize 
time  evolution  of  the  fluorophore  concentration  by  a  slow-time  parameter  f;  and  photon  propagation  during  one 
instance  of  a  tomographic  data  collection  process  by  a  fast-time  parameter  t' .  Note  that  t'  is  in  the  order  of  the 
speed  of  light  whereas  t  is  in  the  order  of  seconds.  Thus,  we  assume  that  the  absorption  and  scattering  coefficients 
of  tissue  are  constant  during  one  instance  of  the  tomographic  data  collection  period  but  vary  with  the  slow-time 
variable  t.  As  a  result,  frequency  domain  couple  diffusion  equation  is  adequate  to  model  light  propagation  during 
the  dynamic  data  collection  process.  Figure  2  illustrates  a  typical  time  evolution  of  the  fluorophore  concentration 
with  respect  to  slow-time  variable  t. 
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In  the  following  subsections,  for  notational  brevity,  we  drop  the  slow-time  dependence  in  our  equations  and  set 
a(r,t )  =  a(r),  paxfir,t)  =  ^(r),  fx,m(r,u,t)  =  fXtm(r,io),  etc.  Note  that  BX)m  is  a  shorthand  notation  for 
the  quantity  B  at  either  excitation  or  emission  wavelengths. 

4)  Model  for  Light  Propagation  in  Fluorescing  Medium:  The  light  propagation  at  the  excitation  and  emission 
wavelengths  can  be  modelled  by  the  following  coupled  diffusion  equations:  [42]: 

-V  •  Dx(r)V<j)x(r,u)  +  (p,ax(r)  +  fx{r,oj)  =  S(r,u),  r  G  O  C  l3  (6) 

-V  •  Dm(r)'V(/)m(r,u)  +  (pom(r)  +  — )  <j>m{r,u)  =  uj)-/paxf(r)  1  •  (7) 

V  C  J  1  +  [ujT(r)\z 

where  the  subscripts  x  and  m  denote  the  excitation  and  emission  wavelengths,  respectively.  0x.m(r,  to)  represents 

the  spatially  varying  optical  field  in  the  medium;  ui  denotes  the  modulation  frequency  of  the  source  (the  Fourier 

transform  with  respect  to  the  fast-time  variable  t')\  c  is  the  speed  of  light  inside  the  medium  O;  r(r)  is  the 

fluorophore  lifetime;  7  is  the  fluorophore’s  quantum  efficiency;  f.iax.m(r)  stands  for  the  spatially  varying  absoiption 

coefficient  of  the  medium  at  the  excitation  and  emission  wavelengths,  respectively,  paxf(r )  is  the  fluorophore’s 

absoiption  coefficient;  and  7 paxf{r)  is  the  fluorophore  yield;  S(r,u)  is  the  excitation  source,  Dxrn{r)  is  the 

spatially  varying  diffusion  coefficient  given  by  Dx^m{r)  =  7^ - 7 - — — y,  where  p!sxm(r)  is  the  reduced 

scattering  coefficient. 

The  optical  coefficients  at  the  excitation  and  emission  wavelengths  are  due  to  both  the  endogenous  chromophores 
and  exogenous  fluorophore.  Thus, 

FaX{p)  =  Faxei'f)  T  Pax/(^”)>  (8) 

Fam{f)  —  Pame(^*)  T  Famfix')-  (9) 

Here,  the  subscript  e  denotes  the  endogenous  chromophores,  and  the  subscript  /  denotes  the  exogenous  fluorophore. 
We  choose  Robin-type  boundary  conditions  given  as  [46] 

2  Dx(r)  ^  +  Pfx(r,  u)  =  0,  redfi  (10) 

2  Dm(r)  ^  +  pfm{r,  u)  =  0,  (11) 

where  <90  denotes  the  boundary  of  O,  v  denotes  the  outward  normal  of  the  boundary  <90,  p  is  a  constant  accounting 
for  the  refractive  index  mismatch  between  the  two  regions  separated  by  OO. 
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Let  o(r)  denote  the  concentration  of  fluorophores  at  r  £  Q.  The  relationship  between  a(r)  and  the  fluorophore 
absorption  coefficient,  na{x,m)f  >  is  given  by  [43] 

Va(x,m)f(r)  =  In  10ex, ma{r)  r  €  Cl,  (12) 

where  ex>m  denotes  the  fluorophore  extinction  coefficients  at  the  excitation  and  emission  wavelengths,  respectively. 

5)  Non-linear  Cone entration-to-Measurement  Map:  Let  <l>m(rd,  rs;  to)  denote  the  ratio  of  the  emission  data  to 

the  excitation  data  (normalized  Born  data  [46])  at  the  emission  wavelength  at  the  detector  location  rd  due  to  an 

excitation  source  at  rs.  The  relationship  between  <l>m(rd,  rs;  cu)  and  a(r)  is  given  by 

^  ,  In  106^7  f  ...  1  —  jurir)  ,  .  ■> 

<$>m(rd,rs;u)  =  — - - - r  /  Gm  rdlr;w  a(r  — - — —^^x(r,rs-u;)dr 

<t>x{rd,rs;uj)  Jn  1  +  [wr(r)]2 

=:  J7rd,rs{a)  (13) 


where  (f>x(r,  rs;  uj)  is  the  photon  density  at  location  r  due  to  a  source  at  rs  at  the  excitation  wavelength;  Gm(rd,  r;  uj) 
is  the  Green’s  function  of  (7)  and  (11)  at  rd,  due  to  a  point  source  at  r,  and  Trd.r,  («)  is  the  non-linear  operator 
that  maps  the  total  fluorophore  concentration  a  to  the  normalized  measurements  4>m.  Note  that  the  non-linearity  in 
(13)  is  due  to  the  pamf  dependence  of  Gmird,  t\ uj),  and  p.a(m)X)  dependence  of  (f>x. 

We  assume  that:  i)  The  endogenous  absorption  coefficients  at  the  emission  and  the  excitation  wavelengths  are 
approximately  equal,  i.e.  pame  ~  llaxe-  This  is  a  valid  assumption  for  a  number  of  different  applications  such  as  small 
animal  and  breast  imaging  [47],  ii )  The  diffusion  coefficients  at  both  the  excitation  and  emission  wavelengths  are 
independent  of  the  endogenous  and  exogenous  absoiption  coefficients,  i.e.  Dx  m(r)  ~  1/(3 p!sx  m(r)).  Furthermore, 
the  diffusion  coefficients  are  known  but  can  be  spatially  varying,  iii )  The  lifetime  parameter,  r(r),  r  €  fl,  is  known, 
and  not  necessarily  constant. 

Let  Nd  and  Ns  denote  the  number  of  detectors  and  sources,  respectively.  Let  $  be  the  measurement  vector 
formed  by  concatenating  the  measurements  for  each  source-detector  pair  as  follows: 

$  =  [<f>m{rdl,rSl;u),---  ,$m(rdjM,rSl;w),  ■  ■  ■  ,  <S>m(rdNd,  rSNp,  u)]T .  (14) 


Using  (13)  and  (14),  we  form  the  following  relationship  between  a  and  <f>: 

<&  =:  IF  (a), 


(15) 


where  T  is  an  operator  with  matrix  kernel  in  which  F,:]  =  Trd  r  for  i  =  1,  •  •  •  ,  Nd,  j  =  1,  •  •  •  ,  Ns. 

Note  that  for  notational  brevity,  we  assume  a  single-frequency  measurement  model.  A  multi-frequency  measure¬ 
ment  model  is  a  straightforward  extension  of  the  single-frequency  measurement  model. 
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6)  Linear  Concentration-to-Measurement  Map:  An  approximate  linear  map  between  the  normalized  measure¬ 
ments  and  the  total  fluorophore  concentration  can  be  obtained  based  on  the  assumption  that  the  presence  of 
exogenous  fluorophores  does  not  change  the  optical  coefficients  pax,m  and  Dx,m  [42].  This  assumption  leads 
to  the  following  relationship  between  the  measurements  and  the  fluorophore  concentrations: 


®m{rd,rs-,u)  = 


In  10ex7 


< i>%(rd,rs-,u )  Jn 
=:  Wrd,r3(a), 


Gem{rd,  r;  u)a{r)  *  +  r UJ)d?r 


(16) 


where  G^rd,  r;  to)  is  the  Green’s  function  of  (7)  and  (1 1)  when  =  /rame,  4>%{r,  rs;  oj)  is  the  predicted  optical 
field  at  r  due  to  a  source  located  at  rs  when  pax  =  naXe,  and  Wrd,rAa)  denotes  the  linear  operator  that  maps  the 
normalized  measurement  at  rd  due  to  a  source  at  rs  to  the  total  fluorophore  concentration. 

Forming  a  measurement  vector  as  in  (14),  we  write 


=:  W(a) 


(17) 


where  W  is  the  linear  operator  with  a  matrix  kernel  in  which  W %j  =  Wrd.,r3  for  i  =  1,  •  •  •  ,  Nd  j  =  1,  •  •  •  ,NS. 

Note  that  the  reconstruction  of  the  pharmacokinetic -rate  images  that  will  be  discussed  in  the  subsequent  sections 
is  not  tied  to  any  specific  linearization  method.  Alternatively,  a  different  linear  approximation  to  T  can  be  obtained 
by  computing  its  first-order  Frechet  derivative  to  yield  a  relationship  between  measurements  and  total  fluorophore 
concentration.  Flowever,  we  adopted  the  model  in  (17)  for  image  reconstruction  using  in  vivo  breast  data. 

7)  Pharmacokinetic-rates-to-Measurement  Map:  In  this  section,  we  combine  the  CTM  map  with  the  spatially 
resolved  compartmental  model  to  obtain  a  mapping  between  the  spatially  resolved  pharmacokinetic-rates  and 
sequence  of  boundary  measurements.  We  call  this  composite  map  the  pharmacokinetic-rate-to-measurement  (PTM) 
map. 

Recall  that 

C(r,t)  =  K(an(r))C(r,t)  te[T0,L \],  reQ.  (18) 

Let  <F(i)  denote  the  measurement  vector  (16)  and  a(r,  t)  denote  the  flourophore  concentration  at  slow-time 
parameter  t.  Combining  the  CTM  map  (15)  with  the  compartmental  model  (2),  we  obtain  the  following  non-linear 
relationship: 

$(f)  =  ^(V(an(r))C(r,f)),  i€[T0,Ti],  r  e  Q.  (19) 

For  the  linear  CTM  map,  (19)  becomes 


*(t)  =  W(V(an(r))C(r,t))  tG[T0,Ti],  reft. 


(20) 
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The  equations  (18)  &  (19)  and  (18)  &  (20)  constitute  the  PTM  map.  We  next  discretize  these  equations  and 
incorporate  dynamic  model  uncertainties  and  measurement  noise  to  the  PTM  map. 

We  use  first-order  Lagrange  basis  to  discretize  the  domain  of  interest  (L  Let  rj,  j  =  1. ....  N  be  the  discrete 
points  representing  the  spatial  location  of  the  voxels  in  Cl.  Let  C(rv  t)  represent  the  concentration  vector  at  time 
t  in  different  compartments,  and  otn  (r; )  represent  the  pharmacokinetic-rates  and  volume  fractions  at  the  jth  voxel 
centered  at  r? ,  j  =  1, ...,  N.  Assuming  that  the  dynamic  measurements  are  collected  at  time  instances,  t  =  fcA, 
k  =  1. ....  K,  where  A  is  the  sampling  period,  we  define  C(r?,  k)  =  C(rv  A: A),  and  express  the  discrete  spatially 
resolved  compartmental  model  as  follows: 


C(rj,k  +  l)  =  K(On(rj))C{rj,k)  +  £{rj,k),  k  =  l,...K,  j  =  (21) 

where  £(r_j,  k)  is  a  zero-mean  Gaussian  process  with  E[£(rj,  k\  )£(r,,  k2)\  =  S(rj  —  ri)5(k\  —  k2)Qn,  representing 
the  dynamic  model  uncertainty;  K(9n(rj))  :=  )  >A  js  qlc  discrete-time  system  matrix  as  described  in  [48] 

and,  0n{rj)  represents  the  discrete-time  parameter  vector  for  the  pharmacokinetic-rates  and  volume  fractions.  For 
a  detailed  discussion  of  the  discretization  procedure  and  an  explicit  relationship  between  the  parameters  On(rj) 
and  OLn(rj),  see  [37],  [48]. 

Let  \P( k )  =  4?(A;A).  Replacing  an(rj)  with  dn(rj),  (19)  and  (20)  are  discretized  as: 


*(k)  =  F(V(0n(rj))C(rj,  k))  +  r/(k),  (22) 

*(k)  =  W\ (6n(r j))C(r j ,  k)  +  r}(k),  k=l,...K,  j  =  l,...N,  (23) 


where  F  and  W  are  the  resulting  operators  when  an(rj)  is  replaced  with  6n(rj)  in  (13)  and  (16),  respectively; 
r}(k)  is  a  zero-mean  Gaussian  process  with  E(ij(k\ ) r/ ( A:-2 ) ]  =  5(k\  —  ) R  representing  the  measurement  noise, 

and  \{Qn(rj))  is  the  vector  of  discrete  volume  fractions  at  rj. 

For  the  two-compartment  ICG  pharmacokinetic  model  combined  with  the  linear  CTM  map,  the  explicit  form  of 
(21)  and  (23)  are  given  by 


Ce(r j,  k  +  1) 
Cp(r j ,  k  +  1) 


Tll(rj) 

1-Mrj) 

1 

■^-3 

1 _ 

+ 

1 

as 

'■T 

_ i 

r2l(ri) 

T22 (rj)  _ 

C p{rj )  k) 

£p(rj >  k) 

(24) 


’F(fc)  =  W 


ve(rj)  vp(r  j) 


Ce(rj,k) 

Cp(rj,k) 


+  V(k),  k  =  l,...,K 


(25) 


where  C (rj,  k)  =  [ Ce(rj,k )  Cp(rj,  k)]T,  £(rj,  k)  =  [£e{rj,  k)  £p(rj,  k)]T,  Y(02(rj))  =  [ve(rj)  vp(rj )],  02(rj) 
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[rn (rj)  T12(rj)  721  (rj)  t22 (rj)  ve(rj)  vp(rj)]T,  and  K (On{rj))  is  the  2-by-2  system  matrix  in  (24),  for 
j  =  1, N  and  k  =  1, K. 

8)  Reconstruction  of  Pharmacokinetic  Rate  and  Concentration  Images  from  the  Boundary  Measurements:  The 
forward  model  in  (21)  and  (22)/(23)  forms  a  state-space  model,  (21)  being  the  state  equation  and  (22)/(23)  being 
the  measurement  equation.  In  this  section,  we  discuss  the  estimation  of  the  system  parameters,  On(rj),  and  the 
states  C d(rj,  k)  for  j  =  1, . . . ,  N,  from  the  measurements  (  /,:),  k  =  1, . . . ,  K  using  the  extended  Kalman  filtering 
(EKF)  framework. 

Note  that  both  the  fluorophore  concentrations  in  different  compartments,  C(r?,  k),  and  the  system  parameters, 
6n(rj),  are  unknown.  In  this  case,  we  estimate  both  the  states  and  system  parameters  from  measurements  within 
the  EKF  framework.  To  do  so,  we  regard  the  state  equation  in  (21)  as  a  non-linear  equation  in  which  the  system 
parameters  and  states  are  combined  to  form  the  new  states  of  the  non-linear  equation.  We  then  iteratively  linearize 
the  non-linear  state  equation  and  solve  for  the  new  unknown  states  using  the  EKF  framework.  This  approach 
requires  use  of  temporal  prior  models  on  Qn(rf).  We  describe  one  such  model  in  the  following  subsection. 

9)  A  priori  Model  for  Pharmacokinetic-rates  and  Volume  Fractions:  To  impose  a  temporal  prior  model  on  On(rj), 
we  extend  our  notation  to  6n(rj,  k),  k  =  1, . . . ,  K.  Note  that  9n(rj,  k )  is  a  vector  containing  pharmacokinetic-rates 
and  volume  fractions  at  location  r3  e  Q  and  time  kA.  For  each  element,  9n(rj,k),  of  the  vector  0n(r1:  />’),  we 
impose  the  following  dynamic  model: 

On(rj,  k+  1)  =  9n{rj,  k)  +  <fi (rj,k)  (26) 

where  (rj,  k )  is  a  zero-mean  Gaussian  process  with  ki)fi(r2,  k2 )]  =  5{r\  —  r2)5(ki  —  k2)z±,  z  1  >  0. 

Note  that  On(rj )  is  modeled  as  a  time-independent  parameter,  and  the  model  in  (26)  relates  9n(rj,k  +  1)  and 
9n{rj,k )  with  an  all-pass  filter.  If,  on  the  other  hand,  6n{rj)  is  time-dependent,  a  different  filter  can  be  chosen 
based  on  a  priori  physiological  information  and/or  robustness  considerations. 

In  addition  to  the  temporal  prior,  we  impose  a  spatial  smoothing  prior  on  0n(rj,k)  to  improve  the  robustness 
of  the  reconstruction  with  respect  to  measurement  noise  and  to  incorporate  a  priori  physiological  information  into 
image  reconstruction.  This  model  is  given  as 

M 

9n{p  j-,  k)  =  ^  Pi9n(rjnk)  +f2{rj,k),  j  =  l,...,JV  (27) 

where  ji,  l  =  1, . . . ,  M  ( ji  f  j )  are  the  indices  of  the  voxels  in  the  neighborhood  of  the  jth  voxel;  /)/,/  =  1 .... .  M 
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are  the  spatial  weighting  coefficients,  which  may  be  different  for  each  pharmacokinetic-rate  or  volume  fraction 
image;  and  ^(r^,  k),  is  a  zero-mean  Gaussian  process  with  -E[si(ri,  fet)^(r2,  ^2)]  =  <5(ri  —  r2)S(k\  —  k2)z2, 
Z2  >  0. 

The  weighting  coefficients  / %  may  be  spatially  varying  and  can  be  chosen  based  on  a  variety  of  physiological 
information,  i.e.,  tumor  location,  size  or  shape.  In  our  numerical  simulations  and  in  vivo  data  processing,  we  assumed 
that  no  such  specific  information  about  the  tumor  is  available  and  used  equal  weights,  i.e.,  (3i  =  1/M.  This  choice 
imposes  an  isotropic  smoothing  on  the  6n(rj,k)  estimates. 

Inserting  the  right-hand  side  of  (27)  for  9n(rj,k )  in  (26),  we  obtain  the  following  spatio-temporal  model  for 
each  entry  dn(r?,  k)  of  the  vector  0,,(r:j.  A:): 


M 

9n(rj,k+ 1)  =  ^2  Pi0n(rjnk)  +  g(rj,k),  k  =  l,...K,  j  =  l,...N,  (28) 

where  g(rj,  k)  is  a  zero-mean  Gaussian  process  with  E[g(r  1,  k\)g{r 2,  fo)]  =  S(r  1  —  r2)6(k\  —  k2)z,  z  >  0,  and  z 
is  a  function  of  z\  and  22-  Figure  ??  illustrates  the  resulting  neighborhood  system  for  2D  images  for  M  =  4. 

Note  that  it  is  possible  to  develop  an  alternative  spatio-temporal  neighborhood  system  taking  into  account  the 
4D  nature  of  the  parameter  space. 

To  simplify  our  notation,  we  express  (28)  in  vector  notation  for  all  entries  of  the  vector  On(rj,k )  as  follows: 


Gn(rj,k+  1)  =  0j(0n(rjl,k),...,0n(rjM,k))  +  <;(rj,k),  k  =  j  =  1 


(29) 


where  /3  -  is  a  vector-valued  linear  function  of  9n{rn ,  /c) ,  l  =  1  ...,M  as  defined  in  (28)  and  q(rj,k)  is  formed  by 
concatenating  the  ?( rj,k )  into  a  column  vector.  It  is  a  zero-mean  Gaussian  process  with  E[s(n,  fci)c(j’2,  fo)]  = 
5(r\  —  r2)5(ki  —  k2)1‘i. 

10)  Estimation  of  Pharmacokinetic-Rate  Images  by  Extended  Kalman  Filtering:  Our  objective  is  to  estimate  the 
fluorophore  concentration  images  in  different  compartments,  and  pharmacokinetic-rate  images.  To  do  so,  we  first 
concatenate  the  concentration  vectors  C (rj,  k)  and  the  parameter  vectors  On(r:j,  k)  for  all  voxels,  j  =  1, ...,  N  and 
form  the  following  vectors: 


C  (A;) 


C  T(ri,k) 


T 


C  T(rN,k)  , 


0n(k) 


On{ri,k)  ...  0*(rN,k) 


T 
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rj>  k+1 

rj2’k 

Fig.  3.  An  illustration  of  the  spatio-temporal  neighborhood  system  for  2-D  pharmacokinetic-rate  and  volume  fraction  images  for  M=A. 
Based  on  the  model  in  (28)  to  (30)  the  neighborhood  of  the  voxel  at  rJ:k+i  (gray  square)  is  given  by  the  voxels  at  rjltk,  rj2:k,  rj3,k,  and 
rjitk  denoted  by  gray-dashed  squares. 


Next,  we  concatenate  the  vectors  C (k)  and  0n(k)  and  form  the  new  state-space  model  based  on  (21)  and  (29). 

(30) 


where 


C(AH-l) 

_ 

K(0n(fc))C(A;) 

+ 

m 

@n(k  +  1) 

(3(en(k)) 

.  ^ . 

'S?(k)  = 


<&(k)  = 


F(V(0n(fc))C(fc))  0 


lnWexyV\(6n(k))C(k)  0 


i 

i 

9n(k) 

+  v(k) 


1 

1 

0n{k) 

+  V(k),  k  =  l,...,K 


K(0n(k))C(k)  = 
P(0n(k))  = 
F(V(0„(A:))C(fc))  = 

m  = 

s(k)  = 


K(9n(r1,k))C(r1,k)  ...  K(9n(rN,k))C(rN,k) 

T 

k), ...)  ...  pN(en(rN,k),...) 


F(V(0„(n  ,k))C(n,k))  ...  F(\(0n(rN,k))C(rN,k )) 

1 T 

£(ri,k)  ■■■  £(rN,k) 

T 


<i{ri,k)  ...  s(rN,k) 


(31) 


(32) 


(33) 


where  £(k)  and  q(k)  are  zero-mean  Gaussian  processes  with  covariance  matrices  Q  and  Z,  respectively.  rj(k)  is 
the  zero-mean  Gaussian  process  with  covariance  matrix  R  as  defined  before. 


16 


Note  that  although  (3  is  linear  in  6n(rj,k),  j  =  1, ...  Ar,  (30)  is  non-linear  in  new  states.  Furthermore,  the  state 
equation  (30)  is  not  block  diagonal. 

We  next  utilize  the  EKF  machinery  to  estimate  the  fiuorophore  concentration  images  in  different  compartments 
and  pharmacokinetic -rates  based  on  (30)  to  (33).  Table  I  tabulates  the  steps  of  the  EKF  algorithm.  The  terms  in 
Table  I  are  defined  as  follows:  C(k\k  —  1)  and  6 n ( k \ k  —  1)  are  the  concentration  and  parameter  estimates  at  time 
k  given  all  the  measurements  up  to  time  k  —  1,  respectively.  Similarly,  C (k)  and  6n(k)  are  the  concentration  and 
parameter  estimate  updates  at  time  k,  respectively.  Pk,k-i  denotes  the  error  covariance  propagation  at  time  k  given 
all  the  measurements  up  to  time  k  —  1;  P is  the  error  covariance  update  at  time  k.  H/  is  the  recursive  Kalman 
gain  matrix  at  time  k  and  I  is  the  identity  matrix.  Jk-i  is  the  Jacobian  matrix  due  to  iterative  linearization  of  the 
(30)  around  C (k  —  1)  and  0n  [k  —  1)-  A fc| ,  is  the  matrix  formed  by  the  discretized  Frechet  derivatives  of  T  with 
respect  to  C  and  0n  at  the  updates  C(k\k  —  1)  and  0n(k\k  —  1).  In  Table  I,  the  EKF  algorithm  is  presented  for  the 
non-linear  case.  For  the  linear  case,  the  non-linear  operator  F  is  simply  replaced  by  the  linear  operator  W. 

The  first-order  Frechet  derivative  of  kF  (or  kFrd.rj  with  respect  to  the  fiuorophore  concentration  a  at  the  EKF 
total  concentration  estimate  at  the  (k\k  —  l)lh  step  is  given  by  [44]— [46] 


9^Frd,r„{^a)  ~ 


in  10eT 


®kJk  1(rd,r , 


—  ([  Gt~\rd,  r;  rs]  co)  1  ^^Sa^dr 

o)\Jn  l  +  [uT(r)y 


-  [  G^-Hrd,  r;  ^k-\r,  rs-u)  1  ^ Sa(r)dr 

Jn  1  +  [wr(r) \z  ex 

~  Gk^~l(rd,  r- cj)^|fc“1(r,  rs;w)  Sa{r)dr^j  (34) 

where  we  define  1(r,rs;o;)  as  the  solution  of  (45)  and  (10),  4>mfc  1(r,rs;tu)  as  the  solution  of  (7)  and  (11), 
Grn'  1  (r\i,  r;  oj)  as  the  Green’s  function  of  (7)  and  (1 1),  and  Gmx  1{r(i, r-  w)  as  the  solution  of  (45)  and  (10)  where 
S(r,<jj )  in  (45)  is  replaced  by  "/fiax/Gm'  1  (rj.  r;  lo)  given  the  EKF  estimates  of  the  fiuorophore  concentrations  at 
different  compartments  at  the  (k\k  —  l)th  step  [44].  In  (34),  the  first  integral  results  from  the  right-hand  side  of  (7), 
while  the  second  and  third  integrals  originate  from  the  dependence  of  and  //„.K,  respectively,  on  the  unknown 
fiuorophore  absoiption  coefficient.  We  note  that  the  kernels  of  the  second  and  thirds  integrals  are  much  smaller 
than  the  kernel  of  the  first  integral.  Therefore,  the  first  integral  in  (34)  dominates  and  the  rest  can  be  neglected.  As 
a  result,  The  first-order  Frechet  derivative  of  T  (or  Trd.r,)  with  respect  to  the  fiuorophore  concentration  a  at  the 
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TABLE  I 

EKF  ALGORITHM  FOR  THE  SIMULTANEOUS  ESTIMATION  OF  FLUOROPHORE  CONCENTRATION,  PHARMACOKINETIC-RATE  AND  VOLUME 

FRACTION  IMAGES. 


System  Model 

C(fc  +  1) 

K(0„(fc))C(fc) 

m 

= 

+ 

0n{k  +  1) 

P  (ffn(k)) 

4*0 

Measurement  Model 


Initial  Conditions 


State  Estimate  Propagation 


'E(fc)  = 

C(l) 

0n(l) 

C(k\k  —  1) 
6„{k\k  -  1) 


F(V(0n(fc))C(fc))  0 


E(C(1)) 

0n(  1) 


1 

0n(k) 

Pl,l  = 


+  T)(k) 

Var{  C(l))  0 

0  Z 


K(0n(fc  —  l))C(fc  —  1) 
P(On(k-  1)) 


Error  Covariance  Propagation  Pt,n  »  Ji-iPt-i,fc-iJ[_i  + 


State  Estimate  Update 


Q  0 
0  z 


C(k) 

C{k\k-1) 

0n{k) 

P(0n(k\k-1)) 

Error  Covariance  Update 


+H *[¥(fc)  -  F(V(e„(fc[fc  -  l))C(fc|fc  -  1))] 
Pfc.fe  =  [I  —  HfeAfeifc-rJPfe,*;-! 


Kalman  Gain 


Hfc  —  Pfc,*;-i  Ak,k_1  [Afcifc-iP^/c-i  Afeifc_1  +  R] 


Definitions 


Jfc-i  — 


K (0n(k  -  1))  ^(K (6n(k  -  1))C (k  -  1)) 

&P(On{k-  1)  ^  /3(fl„(fc-l)) 


Non-linear  Case 


£F(V(0n(fc[fc  -  l))C(k\k  -  1))  J-F(V(04fc|fc  -  l))C(fc|fc  -  1)) 


Linear  Case 


lnlOexW^(\(en(k\k-  l))C(fc|fc  —  1))  lnlOea;Wj-(V(0„(fc|fc-  l))C(fc|fc-  1)) 
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EKF  total  concentration  estimate  at  the  ( k \ k  —  l)"'  step  can  be  approximated  by 


d^rd,rs{^a)  ~ 


in  107e2 


®xk  1{rd,rs-u)  Jn 


Gk^-\rd,r-u)^k-\r,rs-u:) 


1  -  jwr(r) 


5a{r)dr: 


(35) 


1  +  [wr(r)]2 

Based  on  the  assumptions  i)-iii)  in  CTM  map  subsection,  the  terms  in  the  kernel  of  (35)  can  be  computed  as 
follows:  Based  on  the  assumption  i),  fj,nxxi  ~  Home-  Furthermore,  fiamf  and  jinx]'  are  linearly  dependent  by  (12). 
Thus,  given  the  naxe  reconstruction  based  on  (45)  and  (10)  at  the  kth  step,  and  namj  obtained  via  EKF  estimates 
(fluorophore  concentrations  at  different  compartments  and  volume  fractions)  at  the  (k\k—l)th  step,  1  (rd,  r;  to) 
can  be  computed  using  (7)  and  (11).  Similarly,  given  the  fiax,e  reconstruction  based  on  (45)  and  (10)  at  the  kth 
step,  and  fiaxf  estimate  obtained  via  the  EKF  estimates  of  fluorophore  concentrations  at  different  compartments 
and  volume  fractions  at  the  (A:|  A:  —  l)th  step,  <&}xk  1  (rd,rs;  lo)  can  be  computed  as  the  solution  of  (45)  and  (10). 
Finally,  we  assume  that  is  known. 

Note  that  the  calculation  of  the  Frechet  derivative  can  be  simplified  under  some  additional  assumptions.  If  naxe 
does  not  vary  with  the  slow-time  variable  k,  and  that  fiarnf  <C  /iame,  then  Grn(rd.  r:  uj)  can  be  computed  with 
respect  to  fiaxe  and  remains  invariant  with  respect  to  the  slow-time  variable  k. 

Using  the  chain  rule,  the  first-order  Frechet  derivative  of  J7rd,rs  with  respect  to  each  element  of  C  is  given  by 

lnl07€x  f  (7^|fe-i(rdir;a;)$fc|fe-i(r;rs;a;)  1  ^T}r}  SCiVi(r)dr,  (36) 


dFwtfCi)  =  - 

(rd,rs;u)  Jn 


1  +  [wr(r)]2 

where  Vi,i  =  1,2,  ...,n,  denotes  the  volume  fractions. 

Similarly,  the  first-order  Frechet  derivative  of  J7rd,rs  with  respect  to  each  element  of  an  is  given  by 


,r  \  In  1076a 

oSrd,r3{°an)  —  TfcjfcTl - 


Gk [k  1{rd,r-uj)^k  1(r,rs-u)^—^^Svi(r)Ci(r)dr.  (37) 


‘h.r"  (rd,rs;u)  dn 

Note  that  an  contains  both  pharmacokinetic-rate  and  volume  fraction  parameters.  However,  the  Frechet  derivative 
of  T  with  respect  to  pharmacokinetic -rates  is  zero  since  T  depends  only  on  the  volume  fractions. 

In  our  numerical  simulations,  we  used  finite  elements  with  piece-wise  linear  first-order  Fagrange  polynomials  to 
discretize  the  domain  Q  and  thus  to  discretize  the  Frechet  derivatives  given  in  (36)  and  (37). 

The  EKF  algorithms  for  the  estimation  of  concentration  and  parameter  images  are  different  for  the  linear  and 
non-linear  measurement  models  given  in  (13)  and  (16).  When  the  linearized  measurement  model  (16)  is  employed, 
the  operator  W  is  assumed  to  be  constant  throughout  the  dynamic  update  of  the  fluorophore  concentration.  This 
leads  to  the  inherent  assumption  that  the  dynamic  changes  in  the  fluorophore  concentration  can  be  modeled  as  a 
perturbation  on  the  endogenous  chromophore  concentrations.  For  the  non-linear  measurement  model  (13),  on  the 
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other  hand,  the  Frechet  derivative  is  updated  at  every  iteration  of  the  EKF  algorithm  based  on  the  ( k  \  k  —  1 ) th 
update  of  the  concentration  and  volume  fraction  estimates.  While  computationally  more  intense,  the  EKF  algorithm 
based  on  the  non-linear  measurement  model  eliminates  the  limiting  assumption  of  restricting  the  dynamic  changes 
to  perturbations  from  a  constant  endogenous  chromophore  concentration. 

11)  Convergence  and  Initialization  of  EKF:  The  convergence  properties  of  the  EKF  has  been  studied  in  the 
literature  [49],  [50],  [53]-[55].  In  general  for  the  joint  estimation  of  parameters  and  states  the  estimates  may  be 
biased  or  divergent. 

The  main  cause  of  divergence  in  EKF  is  due  to  the  fact  that  there  is  no  coupling  term  between  H/  and  Qn 
[50],  [53].  This  lack  of  coupling  between  H/  and  6n  may  lead  to  divergence  of  the  estimates.  To  overcome  this 
issue,  for  improved  asymptomatic  convergence  properties,  the  EKF  algorithm  can  be  modified  using  the  proposed 
method  given  in  [53].  Flowever,  as  stated  in  [50]  and  [53]  this  modification  imposes  a  high  computational  load 
to  the  algorithm.  Flence,  in  practical  applications  with  large  size  matrices,  instead  of  using  the  proposed  method, 
manual  adjustments  of  the  noise  covariances  matrices  Z,  Q,  and  R  can  be  used  to  improve  the  converge  properties 
of  the  EKF  which  is  called  the  “tuning  of  the  filter”  [50],  [53].  In  our  work,  as  proposed  in  [50]  and  [53],  we 
regal'd  the  initial  values  of  these  matrices  (which  are  nothing  but  the  multiples  of  identity  matrix,  i.e.,  oT)  as  tuning 
parameters. 

The  asymptotic  convergence-rate  of  EKF  depends  the  initialization  of  6n  and  C,  and  proper  selection  of  the 
noise  covariance  matrices  Z,  Q,  and  R  [49],  [50],  [53].  The  covariance  matrix  R  is  the  most  important  term  that 
controls  the  convergence-rate  of  the  EKF  [49].  It  has  been  shown  in  [49]  that  with  an  appropriate  choice  of  the 
covariance  matrix  R,  the  asymptotic  convergence-rate  of  the  EKF  may  be  improved  significantly. 

In  this  work,  we  change  the  values  of  covariance  matrices  R,  Q  and  Z  matrices  and  choose  the  values  that  lead 
to  minimum  norm  of  error  covariance  matrix  within  biological  limits. 

Theoretically,  the  state  estimates  can  be  initialized  at  the  expected  value  of  the  ICG  concentrations,  i.e.  ZT[C(1 )] . 
One  approach  to  the  initialization  of  the  parameters  is  to  utilize  the  state-space  presentation  given  in  (21)-(23). 
Since  .E(\I/(1))  =  eWVd(0ri(l))i?[C(l)],  \F(1)  —  eWV(0n(l)).E,[C(l)]  is  a  zero-mean  random  vector.  If  we  express 
the  variance  of  the  measurement  VF ( 1)  in  terms  of  the  variance  of  C(l)  using  the  measurement  model  in  (23),  and 
solve  for  0n,  we  get  the  estimate  6n(  1)  as  the  most  appropriate  value  for  initialization. 

In  this  work,  we  change  the  initial  values  of  concentrations,  pharmacokinetic -rates,  volume  fractions  within 
biological  limits  and  pick  the  values  that  lead  to  minimum  norm  of  error  covariance  matrix. 
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In  depth  discussion  on  the  convergence  properties  of  the  EKF,  and  the  initialization  of  parameters  and  covariance 
matrices  can  be  found  in  [37],  [49],  [50],  [53]-[55]. 

12)  Computational  Complexity  of  the  EKF  based  Reconstruction  Algorithms:  In  this  subsection,  we  derive  the 
computational  complexity  of  the  EKF  based  direct  reconstruction  algorithms  based  on  the  linear  and  non-linear 
models  in  (31)  and  (32)  under  the  assumptions  outlined  in  sections  VI- A  and  VI-B.  We  next  derive  the  computational 
complexity  of  the  EKF  based  voxel-by-voxel  pharmacokinetic-rate  and  concentration  image  reconstruction  algorithm 
that  we  introduced  in  [41]  and  compare  it  with  the  computational  complexity  of  the  algorithms  introduced  in  this 
paper. 

The  computational  complexity  of  one  recursion  of  the  EKF  is  0(2m2h)  +  0(2mh 2)  +  0(m3)  +  0{h 3),  where 
m  denotes  the  dimension  of  the  measurement  vector,  and  li  denotes  the  dimension  of  the  states  [56]. 

For  the  EKF  based  direct  reconstruction  algorithm  based  on  the  linear  measurement  model  (32),  the  number  of 
states  h  is  n(n  +  2 )N,  where  n  is  the  number  of  compartments,  and  N  is  the  number  of  voxels.  Typical  values 
of  n,  K,  N,  and  m  are  tabulated  in  Table  VI.  Thus,  assuming  that  2 m2n(n  +  2)  ~  N2,  2 mn2(n  +  2)2  ~  N 3/2, 
m 3  ~  N 5/2,  and  n3(n+ 2)3  ~  N,  the  computational  complexity  of  direct  reconstruction  algorithm  for  one  recursion 
of  the  EKF  algorithm  is  given  by  0(N3)  +  0(N3)  +  0(N2)  +  0(N4),  which  is  dominated  by  the  O ( A’4 )  term.  For 
the  non-linear  measurement  model  (31),  the  Frechet  derivative  of  J-  has  to  be  computed  at  every  recursion,  which 
has  a  computational  complexity  of  0(NdN 2  +  NSN2)  [61].  Assuming  that  Nf[  and  Ns  <C  N2,  the  computational 
complexity  of  every  recursion  of  the  EKF  algorithm  using  the  non-linear  measurement  model  is  0(N4).  Hence, 
the  EKF  based  direct  reconstruction  algorithms  using  either  the  non-linear  or  linear  measurement  model  have  the 
computational  complexity  of  0(N4). 

TABLE  II 

Possible  range  of  values  of  the  parameters  used  for  complexity  analysis 

Number  of  compartments,  n  2  —  4 
Total  number  of  voxels,  N  576  —  649 
Size  of  measurement  vector,  m  128  —  256 
Total  number  of  steps,  K  50  —  100 

For  the  EKF  based  voxel-by-voxel  construction  algorithm  that  we  introduced  in  [41],  the  number  of  states,  h,  is 
n(n+2)  and  m  =  1.  In  this  algorithm,  the  absorption  coefficient  images  are  reconstructed  prior  to  pharmacokinetic- 
rate  images.  In  general,  the  computational  complexity  of  this  step  is  0{N3K )  for  a  linear  reconstruction  algorithm. 
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where  K  is  the  number  of  slow-time  steps.  Assuming  that,  2n(n  +  2)  ~  AV1/2,  2n2(n  +  2)2  ~  N,  and  n3(n  +  2)3  ~ 
N 3/2,  the  computational  complexity  of  the  voxel-by- voxel  reconstruction  algorithm  for  one  recursion  of  the  EKF  is 
given  by  0(N3K)  +  0(N3^2)  +  0(N2)+0(N)+0(N5^2),  which  is  dominated  by  the  0(N3K)  term.  For  the  non¬ 
linear  measurement  model,  the  computational  complexity  of  reconstructing  the  absorption  images  of  fluorophore  is 
0(N3Kp )  where  p  is  the  number  of  iterations  performed  in  Born-type  non-linear  iterative  reconstruction  schemes. 
Clearly,  if  the  number  of  slow-time  samples  multiplied  with  the  iteration  number,  Kp,  is  of  higher  order  than  the 
total  number  of  voxels,  N,  then,  the  computational  complexity  of  the  direct  reconstruction  algorithm  is  smaller 
than  that  of  the  voxel-by-voxel  algorithm.  However,  for  the  current  dynamic  imaging  systems,  Kp,  is  smaller  than 
N.  On  the  other  hand,  the  EKF  based  direct  reconstruction  algorithms  offer  a  number  advantages  that  may  justify 
such  an  increase  in  computational  requirements. 

The  algorithm  stores  three  covariance  matrices  Q,  R,  and  Z  with  size  nN  x  nN,  mxm,  and  n(n  +  1)N  x  n(n  + 
l)iV,  respectively.  For  each  iteration,  the  algorithm  stores  a  measurement  matrix  of  size  m  x  1,  stores  and  updates 
the  eiTor  covariance  matrix,  P,  of  size  n(n  +  2)N  x  n(n  +  2)N,  Kalman  gain  matrix,  H,  of  size  n(n  +  2) N  x  m, 
A  matrix  of  size  m  x  n(n  +  2 )N,  and  J  matrix  of  size  n(n  +  2 )N  x  n(n  +  2 )N.  The  algorithm  also  stores  all  the 
updates  of  the  concentrations  and  parameters  which  are  of  size  n(n  +  2)N  x  1. 

13)  Numerical  Simulations  and  Pharmacokinetic-rate  Image  Reconstruction  from  in  vivo  Breast  Data:  We  tested 
the  performance  of  our  approach  using  simulated  data,  and  in  vivo  data  acquired  from  three  patients  with  breast 
tumor.  We  first  present  the  numerical  simulations  and  compare  the  performance  of  direct  reconstruction  algorithm 
(for  both  the  linear  and  non-linear  measurement  models)  with  that  of  voxel-by-voxel  reconstruction  algorithm 
presented  in  [41].  Next,  we  present  the  pharmacokinetic-rate  images  reconstructed  from  in  vivo  breast  data. 

14)  Numerical  Simulations:  We  performed  a  simulation  study  using  the  two-compartment  model  for  ICG 
pharmacokinetics  and  the  light  propagation  model.  Using  physiologically  relevant  values  for  the  pharmacokinetic 
rates,  kin,  kout,  keim,  and  volume  fractions,  ve,  vp,  given  in  Table  III,  we  simulated  the  boundary  measurements, 

A:  =  1, . . . .  K.  for  a  tissue-like  2-D  phantom.  The  maximum  transition  rates  of  kin  and  kout  are  simulated 
at  the  center  of  the  image  and  smoothly  decreased  towards  the  boundaries  based  on  the  results  given  in  [38],  [41]. 
The  fluorescence  quantum  efficiency  and  lifetime  of  ICG  are  assumed  to  be  constant  and  set  to  0.016  and  0.56 
ns,  respectively.  The  modulation  frequency  was  set  to  300  MHz.  The  physical  dimension  of  the  2-D  phantom  was 
chosen  6  cm  by  6  cm.  The  image  domain  was  discretized  into  24  by  24  pixels  each  of  size  0.25cm  by  0.25  cm.  As 
a  prior  model,  we  employed  a  four-pixel  neighborhood  model  (3  =  1/4  due  to  rectangular  nature  of  the  geometry. 


22 


24  sources  and  24  detectors  along  the  boundary  of  the  phantom  were  used  to  generate  simulated  data  as  shown  in 
Figure  4.  The  values  of  keim,  ve  and  vp  in  Table  III  correspond  to  the  average  values  of  the  heterogeneities  from 
the  24  by  24  pixel  phantom  images. 

TABLE  III 

Physiological  values  for  numerical  simulations 


Maximum  value  of  kin 

0.037  sec~ 

Maximum  value  of  kout 

0.029  sec- 

kelm 

0.0054  sec 

Ve 

0.3 

Vp 

0.04 

l-Laxe 

0.05  cm-1 

H'SX 

8  cm-1 

6  cm 


Fig.  4.  The  source  detector  configuration  for  the  numerical  phantom.  Rectangular  shapes  present  the  detectors  and  the  triangular  ones 
present  the  sources. 

We  tested  the  EKF  based  direct  reconstruction  algorithm  based  on  both  the  non-linear  (31)  and  linear  (32) 
measurement  models.  In  the  linear  model,  we  computed  the  matrix  W  based  on  the  background  fxax,m,  and  assumed 
that  it  is  constant  throughout  the  dynamic  update  of  the  fluorophore  concentrations.  For  the  non-linear  model,  we 
updated  the  Frechet  derivative  of  F,  at  every  iteration  of  the  EKF  algorithm  based  on  the  (k\k  —  1 ) th  update  of 
the  concentration  and  volume  fraction  estimates  as  defined  in  (35)-(37).  In  the  update  of  the  Frechet  derivatives, 

we  assumed  that  naxe  is  constant,  and  [iame  Hamf-  As  a  result,  we  computed  Gm(r*d,  r;  cu)  once,  and  did  not 

k\k—l/ 

update  at  every  recursion,  but  updated  the  <fv  (r,  rs;  uj)  at  every  recursion. 

Fig.  5a  and  Fig.  6a  show  the  phantom  images  of  the  pharmacokinetic -rates  kin  and  kout .  Fig.  5b  and  Fig.  6b 
display  the  corresponding  kin  and  kout  images  reconstructed  by  the  EKF  based  direct  reconstruction  algorithm  with 
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Fig.  5. 


♦ 


(a)  Original  phantom  image 


(b)  EKF  based  direct  reconstruction  us-  (c)  EKF  based  direct  reconstruction  us-  (d)  EKF  based  voxel-by-voxel  recon¬ 
ing  the  non-linear  measurement  model,  ing  the  linear  measurement  model.  struction. 

Pharmacokinetic-rate  images  of  ki„  for  three  different  reconstruction  algorithms. 


(a)  Original  phantom  image 


(b)  EKF  based  direct  reconstruction  using  (c)  EKF  based  direct  reconstruction  using  (d)  EKF  based  voxel-by-voxel  recon- 
the  non-linear  measurement  model.  the  linear  measurement  model.  struction. 


Fig.  6.  Pharmacokinetic-rate  images  of  kout  for  three  different  reconstruction  algorithms. 
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the  non-linear  measurement  model.  Fig.  5c  and  Fig.  6c  display  the  corresponding  kin  and  kout  images  obtained 
by  the  direct  reconstruction  with  the  linear  measurement  model.  Figures  5d  and  6d  display  the  corresponding  kin 
and  kout  images  using  the  voxel-by-voxel  algorithm  introduced  in  [41].  We  observe  that  there  is  a  good  agreement 
between  the  true  and  the  estimated  images  in  terms  of  localization  of  the  heterogeneities.  In  all  three  reconstruction 
algorithms,  the  center  of  the  heterogeneity  is  consistent  with  the  ones  in  the  original  phantom  images.  Visual 
inspection  of  reconstructed  images  shows  that  the  EKF  based  direct  reconstruction  algorithm  leads  to  smoother  and 
less  noisy  images  than  that  of  the  voxel-by-voxel  reconstruction. 

To  quantify  the  difference  between  the  reconstructed  and  true  images,  we  used  the  normalized  mean  square  error 


(NMSE): 


NMSE  =  201og10 


|  %  recon  X , 


true  | 


\Xi 


true  | 


Note  that  NMSE  =  —SNR  where  SNR  stands  for  signal-to-noise-ratio.  Table  ??  tabulates  the  NMSE  values  for 
kin  and  kout  images  for  three  different  reconstruction  methods.  For  the  direct  reconstruction  algorithm  with  the 
non-linear  measurement  model,  the  NMSE  for  kin  and  kout  images  are  -19.77  dB  and  -18.49  dB,  respectively.  For 
the  direct  reconstruction  algorithm  with  the  linear  measurement  model,  the  NMSE  for  kin  and  kout  images  are 
-18.45  dB  and  -17.65  dB,  respectively.  Finally,  for  the  voxel-by-voxel  construction  algorithm,  the  error  for  kin  and 
kout  images  are  -16.88  dB  and  -15.90  dB,  respectively. 

Next,  we  studied  the  effect  of  measurement  noise  in  the  performance  of  the  direct  reconstruction  algorithm  and 
compared  it  with  that  of  the  voxel-by-voxel  algorithm.  We  added  zero-mean  white  Gaussian  noise  with  standard 
deviation  equal  to  5%  to  15%  of  the  average  of  the  measurements  with  a  step  size  of  2.5%.  We  generated  15 
realizations  of  the  Gaussian  noise  at  each  level  and  determined  NMSE  based  on  15  reconstructions.  Fig.  7(a)  and 
(b)  show  the  NMSE  versus  the  measurement  noise  for  5  different  noise  levels  for  kin  and  kout.  images,  respectively. 
Clearly,  the  NMSE  in  the  reconstructed  kin  and  kout  images  increases  as  the  measurement  noise  increases.  The 
pharmacokinetic-rate  images  obtained  with  the  direct  reconstruction  algorithm  (both  linear  and  non-linear  cases) 
together  with  the  a  prior  information  results  in  smaller  error  values  as  compared  to  the  voxel-by-voxel  reconstruction 
algorithm.  Moreover,  the  direct  reconstruction  algorithm  using  the  non-linear  measurement  model  resulted  in  smaller 
error  values  as  compared  to  the  direct  reconstruction  algorithm  with  the  linear  measurement  model. 

We  also  studied,  the  effect  of  initialization  of  the  covariance  matrices  R,  Q,  and  Z  in  the  reconstructed  pharmacokinetic- 


rate  images.  We  chose  the  initial  values  of  concentrations,  pharmacokinetic -rates,  volume  fractions  and  R,  Q  and  Z 
matrices  that  lead  to  minimum  norm  of  error  covariance  matrix  within  biological  limits.  Let  R  =  ail,  Q  =  «2l>  and 
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Z  =  C13I,  where  I  is  an  identity  matrix,  and  a.  =  [a\  a 2  0:3].  Table  IV  presents  the  norm  of  the  error  covariance 
matrix  for  kin  and  kout  for  different  values  of  a  using  both  linear  and  non-linear  direct  reconstruction  algorithms. 
The  optimal  a  is  [0.012  0.051  0.0025].  We  observed  that  optimal  a.  results  in  visually  better  quality  images  than 
that  of  other  a  values. 

TABLE  IV 

Norm  of  the  error  covariance  matrix,  a  =  [0.012  0.051  0.0025] 


LINEAR  RECON 

NON-LINEAR  RECON 

a/10 

224.54 

215.93 

a/5 

152.40 

143.36 

a/2 

104.34 

97.34 

(X 

77.18 

63.45 

2a 

107.91 

95.39 

3a 

148.29 

140.70 

In  general,  the  improvements  in  the  direct  reconstructed  images  can  be  attributed  to;  (/)  the  use  of  spatio-temporal 
prior  model  for  the  pharmacokinetic-rate  images  which  leads  to  more  robust  algorithms  and  smoother  images;  (ii) 
efficient  use  of  the  inherent  temporal  correlations  in  NIR  measurements  by  combining  spatial  (photon  propagation) 
and  temporal  (compartmental)  models. 


(a)  (b) 

Fig.  7.  NMSE  vs  measurement  noise  levels  for  the  direct  and  voxel-by-voxel  reconstruction  algorithms  (a)  kin  images,  (b)  kout  images. 


15)  Pharmacokinetic-rate  Images  from  in  vivo  Breast  Data:  We  used  in  vivo  breast  data  acquired  by  a  continuous 
wave  (CW)  NIR  imaging  apparatus  to  reconstruct  the  pharmacokinetic-rate  images  of  ICG.  The  apparatus  has  16 
light  sources  and  16  detectors  located  on  a  circular  holder  at  an  equal  distance  from  each  other  with  22.5  degrees 
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apart.  Sources  and  detectors  were  collocated  and  were  in  the  same  plane.  The  breast  was  arranged  in  a  pendular 
geometry  with  the  source-detector  probes  gently  touching  its  surface.  A  set  of  measurement  for  each  source  was 
collected  at  every  500  ms.  The  total  time  for  the  whole  scan  of  the  breast  including  16  sources  and  16  detectors  was 
8.8  seconds.  The  detectors  used  the  same  positions  as  the  sources  to  collect  the  light  originating  from  one  source 
at  a  time.  Only  the  measurements  from  the  farthest  1 1  detectors  with  high  SNR  were  used  in  image  reconstruction. 
This  resulted  in  approximately  115  viable  measurements  out  of  256  measurements  collected  at  each  time  instant. 
ICG  was  injected  intravenously  by  bolus  with  a  concentration  of  0.25  mg  per  kg  of  body  weight.  Data  acquisition 
started  before  the  injection  of  ICG  and  continued  for  10  minutes. 

Three  patients  with  different  tumor  types  were  included  in  the  study.  First  case,  Case  1,  is  a  fibroadenoma,  which 
corresponds  to  a  mass  estimated  to  be  1—2  cm  in  diameter,  and  located  1  cm  below  the  skin.  Second  case,  Case 
2,  is  an  adenocarcinoma  corresponding  to  a  tumor  estimated  to  be  2—3  cm  in  diameter,  and  located  approximately 
2  cm  below  the  skin.  Third  case,  Case  3,  is  an  invasive  ductal  carcinoma,  which  corresponds  to  a  mass  estimated 
to  be  3—4  cm  in  diameter,  and  located  2  cm  below  the  skin.  Diagnostic  information  was  obtained  by  biopsy  after 
data  acquisition.  A  more  detailed  explanation  of  the  apparatus,  the  data  collection  protocol  and  tumor  information 
can  be  found  in  [40]. 

We  used  a  two-compartment  model  for  the  ICG  pharmacokinetics  as  described  in  (3),  (4)  and  (5).  We  combined 
the  two-compartment  model  with  the  linear  measurement  model  (16).  For  computational  tractability,  we  used  a 
2-D  diffusion  model  for  both  direct  and  voxel-by-voxel  reconstruction.  (See  [68]— [7 1]  for  a  detailed  discussion  of 
errors  resulting  from  using  2-D  diffusion  model  for  3-D  light  propagation  in  breast  tissue.)  We  made  the  following 
simplifying  assumptions:  The  diffusion  coefficient  Dx  is  constant  and  is  equal  to  0.0416  cm.  The  endogenous 
absoiption  coefficient  at  the  excitation  and  emission  wavelength  are  approximately  the  same,  pame  ~  Haxe-  Thus, 
we  determined  G(m  and  <f>%  based  on  (7)  using  the  excitation  measurements  prior  to  ICG  injection.  We  next  set  the 
left  hand  side  of  (16)  to  excitation  measurements  after  the  ICG  injection  and  reconstructed  two-dimensional  ICG 
pharmacokinetic-rate  and  concentration  images  based  on  (30)  and  (32).  The  resulting  measurement  model  is  known 
as  the  differential  diffuse  optical  tomography  model.  A  more  detailed  description  of  the  model  can  be  found  in 
[40],  [41],  [62],  As  a  prior  model,  we  employed  a  six-pixel  neighborhood  model  (3  =  1/6  due  to  circular  nature 
of  the  geometry.  The  initial  values  for  the  pharmacokinetic-rates  and  covariance  matrices  were  regarded  as  tuning 
parameters  and  the  values  that  lead  to  minimum  norm  error  covariance  matrix  were  chosen  as  initial  values. 

The  resulting  ICG  pharmacokinetic -rate  images  are  shown  in  Fig.  8,  9  and  10.  The  images  show  that  there  is  a  good 
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agreement  with  the  location  of  the  heterogeneity  in  the  images  and  the  physical  location  of  the  tumors.  Our  results 
show  that  the  pharmacokinetic -rates  are  higher  around  the  tumor  region  agreeing  with  the  fact  that  permeability 
increases  around  the  tumor  region  due  to  compromised  capillaries  of  tumor  vessels  [63],  [64],  Additionally,  we 
reconstructed  the  ICG  concentration  images  for  plasma  and  the  EES  compartments.  Figures  11-16  show  the  ICG 
concentration  in  plasma  and  the  EES  for  3  different  time  instants  for  Case  1,  2,  and  3,  respectively.  We  observed 
that  the  ICG  concentrations  in  plasma  and  the  EES  compartments  are  higher  around  the  tumors  agreeing  with  the 
hypothesis  that  around  the  tumor  region  ICG  leaks  out  of  damaged  capillaries  of  tumor  vessels. 

Although  the  number  of  available  patient  data  is  limited,  our  results  indicate  that  the  pharmacokinetic -rate  imaging 
may  provide  new  approaches  to  evaluate  and  improve  breast  cancer  diagnosis,  staging,  and  treatment  monitoring. 
Such  approaches  may  include  extraction  of  new  quantitative  features  from  ICG  pharmacokinetic-rate  images,  and 
statistical  analysis  of  spatial  distribution  of  pharmacokinetic -rates. 


(a)  (b) 

Fig.  8.  Case  1:  Direct  reconstructed  pharmacokinetic-rate  images  of  (a)  kin,  (b)  kout. 


(a) 


(b) 


Fig.  9.  Case  2:  Direct  reconstructed  pharmacokinetic-rate  images  of  (a)  kin,  (b)  kout ■ 
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(a)  (b) 

Fig.  10.  Case  3:  Direct  reconstructed  pharmacokinetic-rate  images  of  (a)  kin,  (b)  kout. 


Fig.  11.  Direct  reconstructed  ICG  concentration  images  in  plasma  for  Case  1  for  (a)  246.4th,  (b)  334.4th,  and  (c)  422.4th  seconds. 


(a)  (b)  (c) 

Fig.  12.  Direct  reconstructed  ICG  concentration  images  in  the  EES  for  Case  1  for  (a)  246.4th,  (b)  334.4th,  and  (c)  422.4th  seconds. 


(a) 


(b) 


(c) 


Fig.  13.  Direct  reconstructed  ICG  concentration  images  in  plasma  for  Case  2  for  (a)  228.8th,  (b)  316.8th,  and  (c)  404.8th  seconds. 
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(a)  (b)  (c) 

Fig.  14.  Direct  reconstructed  ICG  concentration  images  in  the  EES  for  Case  2  for  (a)  228.8th,  (b)  316.8th,  and  (c)  404.8th  seconds. 


(a)  (b)  (c) 

Fig.  15.  Direct  reconstructed  ICG  concentration  images  in  the  plasma  for  Case  3  for  (a)  246.4th,  (b)  378.4th,  and  (c)  510.4th  seconds. 


(a) 


(b) 


(c) 


Fig.  16.  Direct  reconstructed  ICG  concentration  images  in  the  EES  for  Case  3  for  (a)  246.4th,  (b)  378.4th,  and  (c)  510.4th  seconds. 
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B.  AIM  3  -  Evaluation  of  NIR  Optical  Features  for  Breast  Cancer  Diagnosis  using  in  vivo  Patient  Data 

The  SOW  with  regard  to  Aim  3  includes  the  following  specific  tasks: 

•  Task  1.  Determine  statistical  variability  of  each  NIR  feature  inside  and  outside  the  suspected  tumor  in  an 
individual  and  evaluate  the  statistical  significance  of  the  measured  difference  with  the  instrumentation  precision. 
12-18tli  month 

•  Task  2.  Design  statistical  classifiers  to  determine  the  ROC  of  each  NIR  feature  for  an  individual.  18-24th 
month 

•  Task  3.  Evaluate  the  ROC,  positive  predicted  value  (PPV)  and  negative  predicted  value  (NPV)  of  various 
combinations  of  the  NIR  features  for  an  individual.  24-27th  month 

•  Task  4.  Investigate  the  significance  of  the  measured  difference  between  malignant  and  benign  tumor  patient 
groups  for  single  and  combined  NIR  features.  27-30th  month 

This  work  presents  the  evaluation  of  a  set  of  optical  features  extracted  from  in  vivo  near-infrared  spectroscopy 
data  obtained  from  116  patients  with  breast  tumors  for  breast  cancer  diagnosis.  The  in  vivo  data  was  collected  from 
44  patients  with  malignant  and  72  patients  with  benign  tumors.  Three  features,  relative  blood  volume  concentration, 
oxygenation  desaturation  and  the  size  of  the  tumor,  are  used  to  differentiate  benign  and  malignant  tumors.  The 
diagnostic  capability  of  these  features  are  evaluated  using  different  classifiers  including  nearest  mean,  neural 
network,  support  vector  machine,  Parzen,  and  normal  density-based  classifiers.  The  area  under  the  receiver  operating 
characteristics  curve  of  the  nearest  mean  classifier  using  the  three  features  yields  the  best  value  of  0.91.  This  result 
suggests  that  relative  blood  volume  concentration,  oxygenation  desaturation  and  size  information  can  differentiate 
malignant  and  benign  breast  tumors  with  a  relatively  high  precision. 

1)  Apparatus:  In  this  study,  a  continuous  wave  (CW)  near  infrared  spectrometer  (NIRS)  was  used  in  collecting 
the  in  vivo  data  [94].  The  apparatus  includes  a  probe.  The  probe  consists  of  one  multi- wavelength  (730nm,  805nm, 
and  850nm)  LED  as  a  light  source  at  the  center  of  the  probe,  and  8  silicon  diodes  as  detectors  arranged  in  a 
circular  geometry  with  4  cm  radius  as  shown  in  Fig.  17.  The  light  intensity  from  the  detectors  was  adjusted  to  be 
approximately  1  volt  and  calibrated  with  a  phantom  with  known  absorption  and  scattering  coefficients  (pa=0.04  to 
0.07  cm-1  and  ps-  8  cm-1). 

2)  Patients  and  Protocol:  The  in  vivo  data  was  collected  at  two  centers,  the  Abramson  Family  Cancer  Research 
Institute,  Department  of  Radiology  of  the  Hospital  of  University  of  Pennsylvania  (HUP),  and  the  Department  of 
Gynecology  of  Leipzig  University  (DGLU).  HUP  provided  24  patients  with  malignant  and  64  patients  with  benign 
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Fig.  17.  The  NIR  probe  with  a  multi-wavelength  LED  and  8  silicon  diodes  as  detectors. 

tumors.  DGLU  provided  20  patients  with  malignant  and  8  patients  with  benign  tumors. 

The  measurements  were  taken  on  the  tumor-free  contralateral  breast  to  include  the  mirror  image  location  of  the 
suspected  tumor.  The  probe  was  then  transferred  to  the  breast  with  suspected  tumor.  The  detectors  giving  the  largest 
changes  with  respect  to  the  contralateral  breast  were  assumed  to  be  related  to  tumor  and  used  for  diagnosis. 

3)  Optical  Features:  In  this  study,  three  features,  namely,  relative  blood  volume  concentration,  ABV,  oxygena¬ 
tion  desaturation,  A Deoxy,  and  the  size  of  the  tumor,  S,  are  used. 

ABV,  and  A  Deoxy  values  are  obtained  based  on  a  lipid  blood  oxygen  model  [94].  Thus,  the  change  in  B  V 
and  Deoxy  are  relative  to  the  contralateral  breast: 

ABV  =  ABVtumor  —  A BVcontra  (38) 

A  Deoxy  =  A  DeoxytUmor  ~  ADeoxycontra  (39) 

where  ABVtumor,  ABVContra  are  relative  blood  volume  concentration  in  the  tumor  breast  with  respect  to  the 
contralateral  breast,  respectively,  and  A Deoxytumor,  ADeoxyCOntra  are  relative  oxygenation  desaturation  in  the 
tumor  breast  with  respect  to  the  contralateral  breast,  respectively. 

The  relative  blood  volume  concentration,  ABV,  and  the  oxygenation  desaturation,  A  Deoxy,  can  be  approximated 
at  two  different  wavelengths  by 

ABV  oc  7AO.D730  +  AOD85o  (40) 

A  Deoxy  oc  /3AOD730  +  AODg^o  (41) 

where  AOD730,  and  AOD8$o  denote  the  relative  changes  in  optical  density  at  730  nm  and  850  nm,  respectively, 
7  =  0.3,  and  f3  =  1.3  are  the  matching  constants,  and  AOD  is  the  differential  optical  density  given  by 

AOD  =  logj 


(42) 
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where  I  is  light  intensity  after  absorption  and  scattering,  and  Jo  is  the  baseline  light  intensity  obtained  from  the 
contralateral  breast. 

To  demonstrate  the  utility  of  equations  (40)  and  (41),  and  obtain  the  matching  constants,  ABV  and  A Deoxy 
were  calculated  using: 

ABV  oc  A [Hb\  +  A [Hb02]  (43) 

A  Deoxy  oc  A{HbC>2 ]  —  A  [Hb\  (44) 

where  A  [Hb\  and  A  [77602  ]  denote  the  relative  change  in  deoxyhemoglobin  (Hb)  and  oxyhemoglobin  (Hb02)  with 
respect  to  the  collateral  breast. 

The  concentrations  of  A [Hb\  and  A [HbOfi^  were  calculated  using: 

AOD  =  eACAL  (45) 

where  OD  is  the  optical  density,  e  is  the  known  extinction  coefficients  of  Hb,  Hb02,  C  is  the  concentration,  L  is 
the  mean  path-length  of  photons.  Here,  e  «  1  cm-1,  and  A L  =  4  cm  for  a  differential  path-length  factor  of  7-8 
[94],  [95], 

4)  Feature  Analysis  and  Tumor  Classification:  In  this  section,  we  present  the  classifiers,  training  techniques, 
statistical  analysis  of  the  dataset,  and  the  malignancy  differentiation  criteria.  Before  we  begin  our  discussion,  note 
that,  in  our  context,  dataset  is  the  set  which  contains  all  malignant  and  benign  cases;  training  and  test  sets  are 
subsets  of  the  dataset. 

Figure  18  shows  the  distribution  of  the  three  features  for  116  patients.  Blue  circles  indicate  malignant  cases  and 
pink  circles  indicates  benign  cases.  Note  that  the  malignant  class  has  a  large  spread  as  compared  to  benign  class. 

5)  Classifiers:  We  evaluate  the  malignancy  differentiation  capability  of  the  individual  features  and  various 
combinations  of  the  these  features  using  5  different  classifiers,  namely,  Parzen  density-based  classifier  (PAR), 
neural  network  classifier  (NEURC),  support  vector  classifier  (SVC),  normal  densities-based  quadratic  classifier 
(NDC),  scaled  nearest  mean  classifier  (NMSC)  [96]— [  100]. 

In  Parzen  density-based  classification  technique,  given  a  kernel  function,  the  probability  distribution  of  the  training 
set  is  approximated  via  a  linear  combination  of  kernels  and  a  test  set  is  assigned  to  the  class  with  maximal  posterior 
probability  [100]. 

The  support  vector  classifiers  are  based  on  support  vector  machines.  In  support  vector  machine  based  classification, 
each  data  point  in  the  dataset  is  represented  by  a  k-dimensional  vector.  Assuming,  each  data  point  belongs  to  only 
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•  Bening 

o  Malignant 


Fig.  18.  Distribution  of  three  features,  A BV,  A Deoxy,  and  S  for  116  patients.  Blue  circles  indicate  malignant  cases  and  pink  circles 
indicates  benign  cases. 

one  of  two  classes,  the  support  vector  classifier  separate  the  dataset  with  a  ”k  minus  1”  dimensional  hyperplane 
with  maximum  separation  between  the  two  classes.  In  other  words,  the  nearest  distance  between  a  data  point  in 
one  hyperplane  and  a  data  point  in  the  other  hyperplane  is  maximized  [96]— [98]. 

The  NDC  is  based  on  computation  of  a  quadratic  discriminant  functional  for  the  classes  in  the  dataset  using 
normal  densities.  Similarly,  NMSC  uses  a  linear  discriminant  functional  for  the  classes  in  the  dataset  assuming 
equal  class  variances  [100]. 

6)  Classifier  Training  Techniques:  The  classifiers  are  trained  by  using  three  different  training  techniques,  namely, 
hold-out,  n-fold  cross  validation,  and  leave-one-out  techniques  [91],  [99],  [100].  In  hold-out  technique  the  original 
dataset  (malignant  and  benign  cases)  is  split  into  training  and  test  sets  randomly.  The  training  set  is  used  for 
generating  the  classification  model  and  the  test  set  is  used  to  test  the  classification  performance  of  the  classifier. 
To  obtain  more  reliable  results  for  hold-out  training  technique,  we  repeated  the  classification  with  different  sub¬ 
samples,  i.e.,  in  each  repetition,  a  certain  proportion  is  randomly  selected  for  training,  the  rest  of  data  is  used  for 
testing.  The  error  rates  on  different  repetitions  are  averaged  to  yield  an  overall  performance  of  the  classifier.  In  the 
hold-out  technique  different  test  sets  overlap  hence  it  is  not  optimal.  In  order  to  prevent  overlapping,  we  used  the 
n-fold  cross-validation  training  technique  with  training  set  size  equal  to  N  —  N/n  and  test  set  size  N/n.  Here  n 
is  the  number  of  subsets  and  N  is  the  total  number  of  the  cases  including  both  malignant  and  being  cases.  In  this 
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technique,  we  split  the  total  number  of  cases  into  n  =  10  subsets  of  equal  size  and  use  each  subset  for  testing  and 
the  remainder  subset  for  training.  The  results  are  then  averaged  to  get  and  overall  classification  performance.  To 
improve  our  results  further,  we  use  a  special  case  of  cross-validation  technique,  the  leave-one-out  technique,  with 
n  =  N.  In  this  technique,  one  sample  is  used  for  testing,  and  the  remanning  N- 1  samples  are  used  for  training  the 
classifier.  More  detailed  discussions  on  these  classifiers  and  training  techniques  can  be  found  in  [91],  [99],  [100]. 

We  evaluated  the  malignancy  differentiation  capability  of  the  following  individual  (A BV,  A Deoxy,  .S’,)  and 
combined  features  {ABV-ADeoxy,  A  BV-  A  Deoxy- S). 

III.  Statistical  Analysis  of  Clinical  Data 

We  evaluated  the  diagnostic  capability  of  different  combinations  of  optical  features  based  on  receiver  operating 
characteristics  (ROC)  methodology  using  different  classifiers  and  training  techniques  [8 1]— [83].  The  ROC  curve 
is  obtained  by  plotting  the  probability  of  false  positive  rate  versus  the  probability  of  detection.  The  evaluation  of 
classification  method  is  done  using  area  under  the  ROC  curve  (AUC). 

First,  we  evaluated  the  classification  performance  of  all  three  features  combined.  Table  V  presents  the  AUC 
values  for  5  different  classifiers  using  three  different  training  techniques.  The  AUC  values  of  5  different  classifiers 
using  leave-one-out  training  technique  range  from  0.8864  to  0.9098  with  NMSC  performing  the  best. 

TABLE  V 

AUC  VALUES  FOR  DIFFERENT  CLASSIFIERS  FOR  ABV-ADeoxy-S 


NMSC 

PAR 

SVC 

NDC 

NEURC 

Leave-one-out 

0.9098±  0.0065 

0.904 1±  0.0048 

0.901 1T0.0057 

0.8984T0.0038 

0.8864T0.0051 

Cross-validation 

0.8905±0.0044 

0.8863T0.0039 

0.8802T0.0066 

0.8737T0.0029 

0.8709T0.0072 

Hold-out 

0.8839T0.0040 

0.8785T0.0058 

0.8743T0.0027 

0.8705T0.0062 

0.8681T0.0066 

TABLE  VI 

AUC  VALUES  FOR  DIFFERENT  CLASSIFIERS  FOR  ABV-ADeoxy 

NMSC 

PAR 

SVC 

NDC 

NEURC 

Leave-one-out 

0.9001T0.0029 

0.8993T0.0043 

0.8902T0.0026 

0.8908T0.0073 

0.8792T0.0048 

Cross-validation 

0.8896T0.0054 

0.8853T0.0070 

0.8830T0.0019 

0.8847T0.0062 

0.8703T0.0071 

Hold-out 

0.8802T0.0048 

0.8761T0.0055 

0.8712T0.0091 

0.8698T0.0047 

0.8622T0.0063 

Next,  we  evaluated  the  performance  of  the  two  features  measured  by  NIR  spectroscopy,  namely  A  BV  and 
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A Deoxy.  Table  VI  presents  the  AUC  values  for  5  different  classifiers.  The  AUC  values  of  5  different  classifiers 


range  from  0.8792  to  0.9001  using  leave-one-out  training  technique.  Figure  19  show  the  distribution  of  features 


A  BV,  and  A  Deoxy  extracted  from  benign  and  malignant  tumors  together  with  5  different  classifiers. 


ABV 


Fig.  19.  5  different  classifiers  and  ABV-ADeoxy  2-D  data  clustering. 

We  next  evaluated  the  individual  classification  performance  of  the  three  features.  Table  VII  presents  the  AUC 
values  for  5  different  classifiers  for  the  feature  ABV.  The  NMSC  has  the  best  performance  in  terms  of  classification 
with  a  AUC  value  of  0.8817.  Table  VIII  presents  the  results  5  different  classifiers  for  the  feature  A  Deoxy.  The 
NMSC  has  the  best  performance  in  terms  of  classification  with  a  AUC  value  of  0.8787.  Table  IX  presents  the 


Fig.  20.  ROC  curves  for  ABV-ADeoxy-S  and  ABV-ADeoxy  using  NMSC  classifier. 
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results  of  5  different  classifiers  for  the  feature  S.  The  NDC  has  the  best  performance  in  terms  of  classification  with 
a  AUC  value  of  0.5612. 

As  it  can  be  seen  from  Table  V,  the  best  performing  feature  set  is  the  combination  of  the  three  features.  However, 
the  combination  set  of  optical  features,  obtained  using  optical  spectroscopy,  can  differentiate  breast  tumors  with  a 
relatively  high  precision  with  a  AUC  value  of  0.9001  (Table  VI).  Similarly,  optical  features,  A BV  and  A Deoxy,  also 
performed  well  with  AUC  values  of  0.8817  and  0.8787,  respectively  (Tables  VII  and  VIII).  We  can  also  conclude 
from  Table  IX  that,  the  tumor  size  alone  is  not  informative  in  differentiating  benign  and  malignant  tumors. 

TABLE  VII 

AUC  VALUES  FOR  DIFFERENT  CLASSIFIERS  FOR  A  BV 


NMSC 

PAR 

SVC 

NDC 

NEURC 

Leave-one-out 

0.8817±0.0036 

0.8802±0.0060 

0.8794±0.0047 

0.8779±0.0032 

0.8513T0.0051 

Cross-validation 

0.8742±0.0058 

0.8719±0.0069 

0.8715±0.0061 

0.8693±0.0070 

0.8461T0.0038 

Hold-out 

0.8703±0.0021 

0.8657±0.0042 

0.8634±0.0018 

0.8617±0.0081 

0.8399T0.0067 

TABLE  VIII 

AUC  VALUES  FOR  DIFFERENT  CLASSIFIERS  FOR  ADeoxy 

NMSC 

PAR 

SVC 

NDC 

NEURC 

Leave-one-out 

0.8787±0.0044 

0.8764±0.0053 

0.8730±0.0069 

0.8711±0.0032 

0.8491T0.0060 

Cross-validation 

0.8703±0.0032 

0.8688±0.0048 

0.8604±0.0076 

0.8598±0.0059 

0.8412±0.0082 

Hold-out 

0.8656±0.0056 

0.8602±0.0041 

0.8567±0.0058 

0.85043=0.0071 

0.8389T0.0046 

Figure  20  presents,  the  ROC  curves  for  all  three  features,  and  the  best  two  features,  namely,  A  Deoxy  and  A  BV 
using  the  scaled  nearest  mean  classifier.  The  observed  area  under  the  ROC  curve  for  ABV-ADeoxy-S,  and  ABV- 
A  Deoxy  are  0.9098,  and  0.9001,  respectively.  Figure  21  presents  the  ROC  curves  for  individual  features  A  BV, 
and  A  Deoxy  using  the  scaled  nearest  mean  classifier.  The  observed  AUC  for  A  BV,  and  A  Deoxy  are  0.8817  and 
0.8787,  respectively. 

Next,  we  checked  whether  the  AUC  values  obtained  for  different  combinations  of  features  are  statistically  different 
or  not.  If  so,  we  can  declare  that  some  combination  of  features  are  more  informative  than  the  others  for  diagnostic 
puiposes.  To  check  the  statistical  difference  in  performance,  we  set  up  a  two  class  hypothesis  testing  problem.  The 
null  hypothesis  corresponds  to  the  case  of  AUCi  =AUC2  while  the  alternative  hypothesis  corresponds  to  (AUCi 
/  AUC2).  The  hypothesis  test  described  above  is  reexpressed  as  follows: 
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TABLE  IX 

AUC  VALUES  FOR  DIFFERENT  CLASSIFIERS  FOR  S 


NMSC 

PAR 

SVC 

NDC 

NEURC 

Leave-one-out 

0.5123T0.0101 

0.5292±0.0088 

0.5291±0.0074 

0.5612±0.0081 

0.5382±0.0092 

Cross-validation 

0.5057±0.0069 

0.5233±0.0084 

0.5204T0.0069 

0.5585T0.0072 

0.5301±0.0113 

Hold-out 

0.5001T0.0078 

0.5184T0.0061 

0.5178T0.0091 

0.5478±0.0046 

0.5267T0.0082 

:  =>  AUCi 

=auc2 

(No  difference 

in  performance) 

H i  :  ==>  AUCi  /  AUC2  (Statistical  difference  in  performance) 

We  assume  that  the  features  are  jointly  Gaussian  distributed  and  computed  the  ^-statistics  as: 

AUC,  -  AUC2 

z  =  /  2  o  =  (46) 

\Jcr{  +  02  —  2rcricr2 

where  r  is  the  correlation  coefficient  of  AUCi  and  AUC2,  erf,  and  a\  are  the  variance  of  AUCi  and  AUC2, 
respectively. 

Next,  we  calculated  the  p-value  based  on  the  ^-statistics  [92],  [93].  We  reject  the  null  hypothesis  Hq  if  the 
p-value  is  smaller  than  0.05  (a  predefined  significance  level),  otherwise  we  accept  the  null  hypothesis. 

We  first  compared  the  AUC  values  of  combined  features  A BV-  A Deoxy-S,  AUC  {ABV- ADeoxy- S),  with 
AUC  ( ABV-ADeoxy )  obtained  by  using  the  NMSC  classifier  trained  with  leave-one-out  technique.  The  p-value  is 
calculated  to  be  0.144.  Using  the  hypothesis  test  described  above,  we  conclude  that  there  is  no  significant  difference 
in  performance  using  the  tumor  size  information  together  with  the  optical  features,  A  BV  and  A  Deoxy.  We  next 
compared  the  AUC  values  of  AUC(ABU)  and  AVC(  ADeoxy)  obtained  by  using  NMSC  classifier.  The  p-value  is 
calculated  to  be  0.361.  Again,  we  found  out  that  there  is  no  significant  difference  between  the  diagnostic  capability 
of  A  BV  and  A  Deoxy.  We  also  compared  the  values  of  A\JC(  ABV-ADeoxy)  with  that  of  AUC(ABV)  and 
AUC(A Deoxy).  The  p-values  are  0.0371,  and  0.0401.  Using  these  values,  based  on  the  hypothesis  testing,  we  can 
conclude  that  combination  of  optical  features  has  better  diagnostic  capability  then  that  of  individual  features  A  BV 
and  ADeoxy.  Finally,  we  compared  the  AUC  values  of  AUC(A/i  U)  and  AUC(A Deoxy)  with  the  AUC  value  of 
AUC(S)  obtained  by  using  NMSC  classifier.  The  p-values  are  0.0103  and  0.0178,  respectively.  Using  the  hypothesis 
test,  there  is  a  significant  difference  in  the  classification  performance  of  the  optical  features,  A  BV  and  ADeoxy, 
when  compared  to  the  tumor  size  information. 

Using  the  same  hypothesis  testing  problem  defined  above,  we  compared  the  performances  of  different  classifiers 
in  terms  of  the  AUC  values  obtained  using  three  features  trained  with  leave-one-out  technique.  The  AUC  value  of 
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NMSC  classifier  is  statistically  different  from  the  AUC  values  of  PAR,  SVC,  NDC,  and  NEURC  classifiers  with 
p-values  of  0.0312,  0.0274,  0.0201,  0.0113,  respectively.  We  also  compared  the  performance  of  PAR  classifier 
with  SVC,  NDC,  and  NEURC  classifiers.  The  AUC  values  of  PAR  and  SVC  classifiers  are  not  statistically 
different  in  terms  of  classification  of  three  features,  with  a  p- value  of  0.182.  However,  the  AUC  value  of  PAR 
classifier  is  statistically  different  from  the  AUC  values  of  NDC  and  NEURC  classifiers  with  p-values  of  0.0234, 
0.0217,  respectively.  Similarly,  the  AUC  value  of  SVC  classifier  is  statistically  different  from  NDC  and  NEURC 
classifiers  with  p-values  of  0.0118,  0.0152,  respectively.  Finally,  we  compared  the  performances  of  NDC  and 
NEURC  classifiers.  The  AUC  value  of  NDC  classifier  is  statistically  different  from  the  AUC  value  of  NEURC 
classifier  with  a  p-value  of  0.0385. 


Fig.  21.  ROC  curves  for  A BV  and  A Deoxy  using  NMSC  Classifier. 
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IV.  Key  Research  Accomplishments 

1)  We  demonstrated,  for  the  first  time  in  the  literature,  the  value  of  spatially  resolved  pharmacokinetic -rates  as 
opposed  to  bulk-rates  using  in  vivo  breast  patient  data.  This  work  was  published  as  a  journal  paper  in  Physics 
in  Medicine  and  Biology  (1).  The  paper  was  downloaded  more  than  250  times  within  the  first  three  months 
after  its  publication  and  was  ranked  among  the  top  10%  of  all  papers  published  by  the  Institute  of  Physics. 

2)  We  developed  methods  of  reconstructing  pharmacokinetic-rate  images  directly  from  NIR  boundary  measure¬ 
ments  and  applied  our  technique  to  in  vivo  breast  cancer  data.  The  results  show  that  the  technique  is  robust 
and  provides  better  signal-to-noise  ratio  images.  We  submitted  our  results  to  IEEE  Transactions  in  Medical 
Imaging. 

3)  We  further  analyze  the  NIR  parameters  collected  from  116  patients  using  NIR  spectroscopy.  This  work  is 
currently  being  prepared  as  a  journal  paper. 

V.  Reportable  Outcomes 

Complete  list  of  outcomes  is  given  below: 

1)  B.  Alacam,  B.  Yazici,  X.  Intes,  B.  Chance,  S.  Nokia,  “Pharmacokinetic -rate  Images  of  Indocyanine  Green 
for  Breast  Tumors  using  Near  Infrared  Optical  Methods,  Physics  in  Medicine  and  Biology,  Vol.  53,  No:  4, 
pp:  837-859,  2008. 

2)  B.  Alacam,  B.  Yazici,  X.  Intes,  B.  Chance,  Direct  Reconstruction  of  Spatially  Resolved  Pharmacokinetic 
Rate  Images  of  Optical  Fluorophores  from  NIR  Measurements,  submitted  to  IEEE  Transactions  in  Medical 
Imaging. 

3)  B.  Alacam,  B.  Yazici.  B.  Chance,  S.  Nioka,  “Evaluation  of  NIR  Optical  Features  for  Breast  Cancer  Diagnosis 
using  in  vivo  Patient  Data,”  to  be  submitted  to  IEEE  Transactions  in  Biomedical  Engineering. 

4)  B.  Alacam,  B.  Yazici.  B.  Chance,  S.  Nioka,  “Characterization  of  Breast  Tumors  with  NIR  Methods  using 
Optical  Indices,  in  Proc.  of  IEEE  EMBS-29th  Anniversary  Conference,  pp.  5186-5189,  Lyon,  France,  August 
2007. 

5)  M.Guven,  B.  Yazici,  E.  Giladi,  “Effect  of  discretization  error  and  adaptive  mesh  generation  for  simultane¬ 
ous  reconstruction  of  scattering  and  absorption  images  in  diffuse  optical  tomography,  Proceedings  of  SPIE 
Photonics  West,  Vol.  6850,  68500W,  January  2008. 
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6)  L.  Reilly-Raska,  M.  Guven,  B.  Yazici,  “Effect  of  discretization  error  and  adaptive  mesh  generation  for 
fluorescence  diffuse  optical  tomography,  Proceedings  of  SPIE  Photonics  West,  Vol.  6850,  68500B,  January 
2008. 

7)  L.  Zhou,  B.  Yazici,  V.  Ntziachristos,  “Fluorescence  Molecular  Tomography  Reconstruction  with  a  priori 
Anatomical  Information,  Proceedings  of  SPIE  Photonics  West,  vol.  6868,  686800,  January  2008. 

8)  M.  Guven,  B.  Yazici,  K.  Kwon,  E.  Giladi,  X.  Intes,  Effect  of  discretization  error  and  adaptive  mesh  generation 
in  diffuse  optical  absorption  imaging:  Part  II,”  Inverse  Problems,  Vol.  23,  pp:  1135-1160,  May  2007. 

9)  M.  Guven,  B.  Yazici,  K.  Kwon,  E.  Giladi,  X.  Intes,  Effect  of  discretization  error  and  adaptive  mesh  generation 
in  diffuse  optical  absorption  imaging:  Part  I,”  Inverse  Problems,  Vol.  23,  pp:  1115-1133,  May  2007. 

10)  M.  Guven,  B.  Yazici,  E.  Giladi,  and  X.  Intes,  “Adaptive  mesh  generation  for  diffuse  optical  tomography,  4th 
IEEE  International  Symposium  on  Biomedical  Imaging:  From  Nano  to  Macro,  pp:  1380-1383,  2007. 

11)  M.  Guven,  B.  Yazici  and  V.  Ntziachristos,  “Fluorescence  diffuse  optical  Image  reconstruction  with  a  priori 
information,  Proceedings  of  SPIE  Photonics  West  2007,  Vol.  6431,  pp:  643107,  2007. 

12)  Murat  Guven,  who  was  partly  supported  by  this  grant,  received  his  PhD  degree  in  December  2007. 

13)  Burak  Alacam,  who  was  partly  supported  by  this  grant  received  his  PhD  degree  in  June  2008. 

VI.  Conclusions 

In  January  2008,  we  published  our  results  on  ICG  pharmacokinetic-rate  imaging  applied  to  three  in  vivo  patient 
data  with  breast  tumors  in  Physics  in  Medicine  and  Biology.  Our  paper  is  downloaded  more  than  250  times  within 
the  next  three  months  of  its  publication.  We  have  received  a  note  from  Institute  of  Physics  that  our  paper  is  ranked 
among  the  top  10%  among  all  journal  papers  published  by  Institute  of  Physics.  The  reviewer  comments  also  pointed 
out  to  the  novelty  of  our  work. 

We  continued  our  work  on  the  reconstruction  of  pharmacokinetic-rate  images  from  NIR  boundary  measurements. 
Instead  of  the  two-step  decoupled  approach  described  in  our  Physics  in  Medicine  and  Biology  publication,  this  time 
we  studied  the  reconstruction  of  pharmacokinetic-rate  images  directly  from  the  optical  boundary  measurements.  This 
approach  couples  the  reconstruction  of  ICG  concentration  images  and  pharmacokinertic-rate  images.  As  a  result,  it 
takes  advantage  of  inherent  temporal  and  spatial  correlations  between  consecutive  optical  boundary  measurements. 
This  results  in  higher  SNR  than  the  pharmacokinetic-rate  images  obtained  from  absoiption  images.  We  assessed 
the  value  of  this  approach  in  in  vivo  data  obtained  from  three  patients  with  breast  tumors  and  showed  that  they 
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result  in  higher  SNR  than  the  images  obtained  with  the  two-step  decoupled  approach.  We  submitted  our  results  to 
IEEE  Transactions  in  Medical  Imaging. 

We  continued  to  analyze  the  116  breast  cancer  patient  data  obtained  using  NIR  spectroscopy  and  refine  our 
imaging  algorithms.  The  in  vivo  data  was  collected  from  44  patients  with  malignant  and  72  patients  with  benign 
tumors.  We  used  three  features,  relative  blood  volume  concentration,  oxygenation  desaturation  and  the  size  of 
the  tumor,  to  differentiate  benign  and  malignant  tumors.  We  evaluated  the  diagnostic  capability  of  these  features 
using  different  classifiers  including  nearest  mean,  neural  network,  support  vector  machine,  Parzen,  and  normal 
density-based  classifiers.  The  area  under  the  receiver  operating  characteristics  curve  of  the  nearest  mean  classifier 
using  the  three  features  yields  the  best  value  of  0.91.  This  result  suggests  that  relative  blood  volume  concentration, 
oxygenation  desaturation  and  size  information  can  differentiate  malignant  and  benign  breast  tumors  with  a  relatively 
high  precision.  This  work  will  be  submitted  as  a  journal  paper  to  IEEE  Transactions  in  Biomedical  Imaging. 

Murat  Guven  and  Burak  Alacam,  who  were  partly  supported  by  this  grant,  received  their  Ph.D.  degrees  in 
December  2007,  and  June  2008,  respectively. 

We  applied  for  a  “No  Cost  Extension”  for  our  grant  to  complete  the  Aim  4  of  our  project.  The  main  impediment 
to  our  progress  has  been  the  retirement  of  our  collaborator,  Prof.  Britton  Chance  of  University  of  Pennsylvania.  His 
retirement  interrupted  the  MR-NIR  patient  data  collection.  The  novel  concurrent  imager  is  now  taken  over  by  Prof. 
Arjun  Yodh  of  University  of  Pennsylvania.  We  reached  an  agreement  with  Prof.  Yodh  to  collaboratively  collect 
data  in  the  next  one  year  and  analyze  it.  We  are  confident  that  we  will  complete  Aim  4  of  our  project  and  submit 
a  final  report  in  the  next  reporting  cycle. 
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Abstract 

In  this  paper,  we  develop  a  method  of  forming  pharmacokinetic -rate  images 
of  indocyanine  green  (ICG)  and  apply  our  method  to  in  vivo  data  obtained 
from  three  patients  with  breast  tumors.  To  form  pharmacokinetic-rate  images, 
we  first  obtain  a  sequence  of  ICG  concentration  images  using  the  differential 
diffuse  optical  tomography  technique.  We  next  employ  a  two-compartment 
model  composed  of  plasma,  and  extracellular-extravascular  space  (EES),  and 
estimate  the  pharmacokinetic  rates  and  concentrations  in  each  compartment 
using  the  extended  Kalman  filtering  framework.  The  pharmacokinetic-rate 
images  of  the  three  patient  show  that  the  rates  from  the  tumor  region  and 
outside  the  tumor  region  are  statistically  different.  Additionally,  the  ICG 
concentrations  in  plasma,  and  the  EES  compartments  are  higher  around  the 
tumor  region  agreeing  with  the  hypothesis  that  around  the  tumor  region  ICG 
may  act  as  a  diffusible  extravascular  flow  in  compromised  capillary  of  cancer 
vessels.  Our  study  indicates  that  the  pharmacokinetic-rate  images  may  provide 
superior  information  than  single  set  of  pharmacokinetic  rates  estimated  from 
the  entire  breast  tissue  for  breast  cancer  diagnosis. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Near-infrared  (NIR)  diffuse  optical  imaging  offers  several  advantages  over  other  imaging 
modalities  (Arridge  1999,  Boppart  etal  2004,  Gu  etal  2004,  Intes  and  Chance  2005,  Mahmood 
et  al  1999,  Sevick-Muraca  et  al  1997,  Yodh  and  Chance  1995).  NIR  techniques  are  minimally 
invasive,  easy  to  use,  relatively  inexpensive  and  can  be  made  portable.  Moreover,  optical 
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techniques,  when  coupled  with  contrast  agents,  have  the  potential  to  provide  molecular/cellular 
level  information,  which  can  improve  cancer  detection,  staging  and  treatment  monitoring 
(Alacam  et  al  2006,  Cuccia  et  al  2003,  Intes  et  al  2003,  Mahmood  et  al  1999,  Sevick-Muraca 
et  al  1997). 

Among  many  commercially  available  optical  contrast  agents,  only  indocyanine  green 
(ICG)  is  approved  for  use  in  humans  by  the  Food  and  Drug  Administration  (ElDeosky  et  al 
1999,  Hansen  etal  1993,  Shinohara  et  al  1996).  ICG  is  a  blood  pooling  agent  that  has  different 
delivery  behavior  between  normal  and  cancer  vasculature.  In  normal  tissue,  ICG  acts  as  a 
blood  flow  indicator  in  tight  capillaries  of  normal  vessels.  However  in  tumors,  ICG  may  act 
as  a  diffusible  (extravascular)  flow  in  leaky  capillary  of  vessels  (Alacam  et  al  2006,  Cuccia 
et  al  2003,  Ntziachristos  et  al  2000,  Vaupel  et  al  1991).  Therefore,  pharmacokinetics  of  ICG 
has  the  potential  to  provide  new  tools  for  tumor  detection,  diagnosis  and  staging. 

One  approach  to  analyze  pharmacokinetics  of  contrast  agents  is  the  compartmental 
modeling  (Anderson  1983,  Jacquez  1972,  Tornoe  2002).  A  number  of  studies  using 
compartmental  modeling  were  reported  to  show  the  feasibility  of  ICG  pharmacokinetics  in 
tumor  characterization  (Alacam  et  al  2006,  Cuccia  et  al  2003,  Intes  et  al  2003).  Cuccia  et  al 
(2003)  presented  a  study  of  the  dynamics  of  ICG  in  an  adenocarcinoma  rat  tumor  model  using 
a  two-compartment  model.  Intes  et  al  (2003)  presented  the  uptake  of  ICG  in  breast  tumors 
using  a  continuous  wave  diffuse  optical  tomography  apparatus  using  a  two-compartment 
model.  We  recently  introduced  the  extended  Kalman  filtering  (EKF)  framework  to  model  and 
estimate  the  ICG  pharmacokinetics  and  tested  three  different  compartmental  models  for  the 
ICG  pharmacokinetics  using  the  in  vivo  NIR  data  collected  from  Fischer  rats  with  cancerous 
tumors  (Alacam  et  al  2006).  Our  study  suggests  that  the  pharmacokinetic  rates  out  of  the 
vasculature  are  higher  in  edematous  tumors  as  compared  to  necrotic  tumors. 

In  all  the  studies  described  above,  the  pharmacokinetic  rates  are  assumed  to  be  constant 
over  a  tissue  volume  that  may  be  as  large  as  the  entire  imaging  domain.  However, 
pharmacokinetic  rates  are  expected  to  be  different  in  healthy  and  tumor  tissue  as  reported 
in  positron  emission  tomography  (PET)  and  magnetic  resonance  imaging  (MRI)  literature. 
It  was  shown  that  the  spatially  resolved  pharmacokinetic-rate  analysis  provides  increased 
sensitivity  and  specificity  for  breast  cancer  diagnosis  (Mussurakis  et  al  1997,  Su  et  al 
2005,  Sun  et  al  2006).  For  example,  Sun  et  al  (2006)  showed  that  FAU  (l-2'-deoxy-2'- 
fluoro-  /3-D-arabinfuranosyl  urasil,  a  PET  contrast  agent)  accumulation  in  tumor  regions  is 
significantly  higher  when  compared  to  normal  breast  tissue  based  on  pharmacokinetic-rate 
images.  Mussurakis  et  al  (1997)  showed  that  the  pharmacokinetics  of  gadolinium-DTPA  (an 
MRI  contrast  agent)  can  be  used  to  differentiate  between  malignant  and  benign  breast  tumors 
with  a  high  accuracy.  It  has  also  been  shown  that  the  spatially  resolved  image  interpretation 
is  superior  to  the  isolated  use  of  quantitative  pharmacokinetic  rates. 

In  the  area  of  diffuse  NIR  spectroscopy  and  imaging,  a  number  of  studies  on  spatially 
resolved  pharmacokinetic  rates  has  been  reported  (Gurfinkel  et  al  2000,  Milstein  et  al  2005). 
Gurfinkel  et  al  (2000)  presented  in  vivo  NIR  reflectance  images  of  ICG  pharmacokinetics  to 
discriminate  canine  adenocarcinoma  (located  at  0.5-1  cm  depth)  from  normal  mammary  tissue. 
These  images  were  generated  by  a  non-tomographic  technique  using  a  CCD  camera  that  is 
suitable  only  for  imaging  tumors  close  to  surface.  Milstein  et  al  (2005)  presented  a  Bayesian 
tomographic  image  reconstruction  method  to  form  pharmacokinetic-rate  images  of  optical 
fluorophores  based  on  fluorescence  diffuse  optical  tomography.  Numerical  simulations  show 
that  the  method  provides  good  contrast.  However,  no  real  data  experiments  were  presented  to 
study  the  diagnostic  value  of  spatially  resolved  pharmacokinetic  rates. 

In  this  paper,  we  present  a  method  of  forming  pharmacokinetic-rate  images  and  report 
spatially  resolved  pharmacokinetic  rates  of  ICG  using  in  vivo  NIR  data  acquired  from  three 
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patients  with  breast  tumors.  To  the  best  of  our  knowledge,  our  work  is  the  first  study 
presenting  the  pharmacokinetic-rate  images  of  an  optical  contrast  agent  using  in  vivo  breast 
data  based  on  tomographic  techniques.  We  first  develop  a  set  of  spatio-temporally  resolved 
ICG  concentration  images  based  on  differential  diffuse  optical  tomography.  We  model  the 
ICG  pharmacokinetics  by  a  two-compartment  model  composed  of  plasma  and  extracellular- 
extravascular  space  (EES)  compartments.  We  then  estimate  the  ICG  pharmacokinetic  rates 
and  the  concentrations  in  different  compartments  based  on  the  EKF  framework  (Alacam 
et  al  2006).  We  show  that  the  pharmacokinetic  rates  from  the  tumor  region  and  outside  the 
tumor  region  are  statistically  different.  We  also  estimate  a  single  set  of  pharmacokinetic  rates 
(bulk  pharmacokinetic  rates)  for  the  entire  breast  tissue.  Our  study  indicates  that  spatially 
resolved  pharmacokinetic  rates  provide  more  consistent  and  superior  diagnostic  information 
as  compared  to  the  bulk  pharmacokinetic  rates. 

The  rest  of  the  paper  is  organized  as  follows.  In  section  2,  we  present  the  reconstruction 
of  ICG  concentration  images.  In  section  3,  we  present  modeling  and  estimation  of  ICG 
pharmacokinetic -rate  images  using  the  EKF  framework.  In  section  4,  we  present  the  spatially 
resolved  ICG  pharmacokinetic-rate  analysis  of  in  vivo  breast  data.  Section  5  summarizes  our 
results  and  conclusion. 

2.  Reconstruction  of  bulk  ICG  concentration  images 

In  our  data  collection  process,  a  sequence  of  boundary  measurements  are  collected  over  a  period 
of  time.  Each  set  of  measurements  are  used  to  form  a  frame  of  the  ICG  concentration  images. 
The  resulting  sequence  of  ICG  concentration  images  are  then  used  to  form  pharmacokinetic- 
rate  images.  To  reconstruct  each  frame  of  the  ICG  concentration  images,  we  follow  a  static 
reconstruction  approach  and  use  differential  diffuse  optical  tomography  (DDOT)  technique 
(Intes  et  al  2003,  Ntziachristos  et  al  1999). 

In  DDOT,  two  sets  of  excitation  measurements  are  collected  corresponding  to  before  and 
after  the  ICG  injection,  and  the  ICG  concentration  is  determined  by  the  perturbation  method 
(Intes  et  al  2003,  Ntziachristos  et  al  1999).  The  photon  propagation  before  and  after  the 
injection  is  modeled  by  the  following  diffusion  equations: 

V  •  ^(rjVit^r,  at)  -  (4(r)  +;(»/c)  <^(r,  a.)  =  0,  r  e  C  R\  (1) 

with  Robin- type  boundary  conditions: 

d<pf(r,  co)  i 

2 Dx(r) — *  +  p<S>f(r,  to)  =  -S(r,co),  r  e  dQ,  (2) 

3v 

where  x  stands  for  the  excitation,  c  is  the  speed  of  light  inside  the  medium  Q;  o>  denotes  the 
modulation  frequency  of  the  source,  fx~x(r)  and  li+ax(r)  are  the  absorption  coefficients  before 
and  after  the  ICG  injection,  Dx  is  the  diffusion  coefficient  which  is  assumed  independent  of 
jx^x,  known  but  not  necessarily  constant,  co)  denotes  optical  field  at  location  r  before 

and  after  the  ICG  injection.  Here,  v  denotes  the  outward  normal  to  the  boundary  3Q  of  Q,  p 
is  a  constant  representing  the  refractive  index  mismatch  between  the  two  regions  separated  by 
3  G  and  S  (r,  co)  is  the  excitation  source  on  the  boundary. 

The  absorption  coefficients  after  the  injection  \±+ax  are  modeled  as  a  sum  of  the  absorption 
coefficient  of  the  medium  before  the  ICG  injection  p,~x  and  the  perturbation  caused  by  the 
ICG  A nax(r)\ 

A n,ax(.r)  =  nlx(r)  -  fx~x(r),  r  e  £2  C  R3.  (3) 

In  the  forward  model,  the  analytical  solutions  of  the  heterogonous  diffusion  equation  given 
in  (1)  is  derived  using  first-order  Rytov  approximation  (Intes  et  al  2003).  The  sample  volume 
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is  divided  into  a  set  of  voxels  and  the  measurements  are  related  to  the  relative  absorption 
coefficients  of  each  voxel  by  a  system  of  linear  equations.  The  shape  of  the  breast  was 
approximated  as  a  cylinder  and  the  Kirchhoff  approximation  (Ripoll  et  al  2001a,  2001b)  for 
diffuse  waves  was  used  to  model  the  interaction  of  light  with  boundaries.  In  order  to  minimize 
optode-tissue  coupling  mismatch  due  to  breathing  motion,  the  forward  model  was  augmented 
with  the  coupling  coefficient  technique  as  described  in  Boas  et  al  (2001). 

Here,  the  Rytov-type  measurements,  which  are  defined  by  the  natural  logarithm  of  the 
ratio  of  the  post-ICG  measurements  to  the  pre-ICG  measurements  were  used  (Ntziachristos 
et  al  1999).  Let  rj;  rs)  denote  the  Rytov-type  measurements  at  location  r,i  due  to 

source  at  rs.  The  linearized  relationship  between  the  differential  absorption  coefficient  and 
measurements  is  given  by  O’Leary  (1996) 


(4) 


where  <t>7(r,  or.  rs)  is  the  photon  density  obtained  at  the  excitation  wavelength  before  ICG 
injection,  Ix(r)  —  cAp,ax(r)/Dx,  and  G~(r,  or,  r,i)  is  the  Green’s  function  of  (1)  for  a  source 
at  rs  before  the  injection. 

We  address  the  inverse  problem  of  recovering  A fiax  from  Rytov  measurements  TA  based 
on  the  forward  model  (4)  using  the  singular-value  decomposition  of  the  Moore-Penrose 
generalized  system.  We  use  a  zeroth-order  Tikhonov  regularization  to  stabilize  the  inversion 
procedure.  The  regularization  parameter  was  determined  by  L-curve  analysis  (Hansen  and 
O’Leary  1993)  using  the  data  obtained  from  a  phantom  study  previously  employed  to  validate 
the  apparatus  (Intes  etal  2003).  The  optimal  regularization  parameter  was  found  to  be  6  x  1 0  4 
and  set  to  be  the  same  for  all  patient  images  and  time  instances.  A  detailed  discussion  of  the 
forward  and  inverse  models  used  for  the  reconstruction  of  differential  absorption  coefficients 
( A fiax )  can  be  found  in  Intes  et  al  (2003). 

To  construct  a  set  of  ICG  concentration  images,  we  use  the  linear  relationship  between 
the  differential  absorption  coefficients  and  ICG  concentrations  (Landsman  et  al  1976), 


A  fia(r)  =  In  10  exm(r), 


(5) 


where  ex  is  the  extinction  coefficient  of  ICG  at  the  wavelength  805  nm,  m  (r)  is  the  bulk  ICG 
concentration  in  the  tissue  and  Afia(r)  is  as  defined  in  (3). 

Note  that  the  method  described  here  is  applicable  for  frequency  domain  case  but  for 
simplicity  we  set  the  frequency  to  zero,  i.e.  a>  —  0. 

3.  Modeling  and  estimation  of  ICG  pharmacokinetics 

3.1.  Two-compartment  model  of  ICG  pharmacokinetics 

Using  the  method  outlined  in  section  2,  we  reconstruct  a  sequence  of  ICG  concentration 
images.  As  an  example,  figures  1-3  show  a  set  of  images  reconstructed  from  in  vivo  breast 
data. 

Our  objective  is  to  model  the  pharmacokinetics  of  ICG  at  each  voxel  of  ICG  concentration 
images  using  compartmental  modeling.  To  do  so,  we  first  extracted  the  time  varying  ICG 
concentration  curves  for  each  voxel  from  the  sequence  of  ICG  concentration  images.  An 
example  of  such  a  curve  is  shown  in  figure  4.  We  next  fit  a  two-compartment  model  to  each 
ICG  concentration  curve  (Alacam  et  al  2006,  Gurfinkel  et  al  2000). 
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Figure  1.  ICG  concentration  images  for  a  set  of  time  instants  for  case  1 . 


Using  the  two-compartment  model  introduced  by  Alacam  et  al  (2006),  ICG  transition 
between  plasma  and  extracellular-extravascular  space  (EES)  can  be  modeled  as  follows: 


'CM 

Cp{t)_ 

^in 

(^in  +  Clm) 


ce(t) 

Cp(t ) 


+  u>(t), 


t  e  [r0,  7)], 


(6) 


where  Cp(t)  and  C,,(t)  represent  the  ICG  concentrations  in  plasma  and  EES  at  time  t. 
respectively.  The  rates  km ,  koat  and  kcim  have  a  unit  of  sec"  1 .  They  are  defined  as  the 
permeability  surface  area  products  given  by  PSy,  where  P  is  the  capillary  permeability 
constant,  S  is  the  capillary  surface  area  and  y  is  the  tissue  density.  hm  and  kout  govern 
the  leakage  into  and  the  drainage  out  of  the  EES.  The  parameter  keim  describes  the  ICG 
elimination  from  the  body  through  kidneys  and  liver.  Here,  u >(t)  is  uncorrelated  zero-mean 
Gaussian  process  with  covariance  matrix  Q  representing  the  model  mismatch. 

The  actual  total  ICG  concentration  in  the  tissue  is  a  linear  combination  of  plasma  and  the 
EES  ICG  concentrations,  and  modeled  as 


lit)  =  [ ve  vp ] 


Ce{t) 

L  Cp(t) 


+  V  (t). 


t  e  [To,  7i], 


(7) 


where  m(t),  Ce(t)  and  Cp(t)  are  defined  in  (5)  and  (6),  vp  and  ve  are  plasma  and  the  EES 
volume  fractions,  respectively,  and  rj(t)  is  uncorrelated  zero-mean  Gaussian  process  with 
covariance  matrix  R,  representing  the  measurement  noise. 
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Figure  2.  ICG  concentration  images  for  a  set  of  time  instants  for  case  2. 


3.2.  Estimation  of  ICG  pharmacokinetics  using  extended  Kalman  filtering 

In  matrix-vector  notation,  (6)  and  (7)  can  be  expressed  as 

C(t)  =  K(a)C(t)  +  w(f),  m{t)  =  V(a)C(f)  +  rj(t),  (8) 

where  C(f)  denotes  the  concentration  vector  with  elements  Ce(t).  and  Cp(t)\  K(a)  is  the 
system  matrix,  V(ct)  is  the  measurement  matrix  as  defined  in  equation  (7)  and  a.  is  the 
parameter  vector  given  by 

tX  =  [koat  k[n  ke  im  t-V;  I'/?  1  ■  (9) 

The  ICG  measurements  in  (8)  are  collected  at  discrete  time  instances,  t  =  kT,  k  = 
0,  1,  . . . ,  where  T  is  the  sampling  period.  Therefore,  the  continuous  model  described  in  (8)  is 
discretized.  We  can  express  the  discrete  compartmental  model  as  follows: 

C  d(k  +  1)  =  Kd(0)Cd(k)  +  u>d(k),  m(k)  =  \d(0)Cd(k)  +  Vd(k),  (10) 

where  K d(G)  —  eK,r>)  is  the  discrete  time  system  matrix;  \ d(0)  =  Y (a)  is  the  discrete 
measurement  matrix;  Lod(k)  and  rjd(k)  are  zero-mean  Gaussian  white  noise  processes  with 
covariances  matrix  Qd  and  variance  R,/,  respectively.  The  vector  6  is  composed  of  parameters 
t ij  which  are  functions  the  pharmacokinetic  rates  and  volume  fractions: 

0  =  [r  11  t!2  r2l  T22  ve  Vp]T . 


(11) 


Y  [cm]  ii  Y  [cm]  ^  Y  [cm] 
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Figure  3.  ICG  concentration  images  for  a  set  of  time  instants  for  case  3. 


Figure  4.  Time  course  of  ICG  concentration  curves  for  a  specific  voxel,  65th,  276th,  188th  voxel 
for  cases  1,  2  and  3,  respectively. 
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Figure  5.  (Left)  Schematic  diagram.  (Right)  The  cut  section  of  the  CW  NIR  imaging  apparatus 
with  16  sources  and  detectors. 


We  first  estimate  t i,  j  =  1,2  and  then  compute  the  pharmacokinetic  rates  £;n,  /cout  and  ke\m 
(Alacam  et  al  2006,  Chen  1999).  The  explicit  form  of  the  discrete  state-space  model  is  given 
as  follows: 


~Ce(k  +  l)' 

"*11 

T12" 

'cm  ' 

Cp(k+  1)_ 

J21 

T22 . 

Cp(k)_ 

+  ud(k) 


m(k )  =  [ve  vp 


'CM 

LCP(k) 


(12) 


+  Vd(k). 


We  estimate  the  parameter  vector  6  and  concentration  vector  Cd  by  using  the  EKF 
framework.  The  EKF  is  a  recursive  modeling  and  estimation  method  with  numerous 
advantages  in  ICG  pharmacokinetic  modeling  (Alacam  et  al  2006).  These  include  effective 
modeling  of  multiple  compartments,  and  multiple  measurement  systems  in  the  presence  of 
measurement  noise  and  uncertainties  in  the  compartmental  model  dynamics,  simultaneous 
estimation  of  model  parameters  and  ICG  concentrations  in  each  compartment,  statistical 
validation  of  estimated  concentrations  and  error  bounds  on  the  model  parameter  estimates,  and 
incorporation  of  available  a  priori  information  about  the  initial  conditions  of  the  permeability 
rates  into  the  estimation  procedure. 

When  both  states  (ICG  concentrations)  and  model  parameters  (pharmacokinetic  rates 
and  volume  fractions)  are  unknown,  a  linear  state-space  model  can  be  regarded  as  a  nonlinear 
model;  the  linear  system  parameters  and  states  combine  to  form  the  new  states  of  the  nonlinear 
model.  This  system  is  then  linearized  and  the  new  unknown  states  are  found  using  the  EKF 
estimator  (Alacam  et  al  2006,  Ljung  1979,  Togneri  and  Deng  2003,  Nelson  and  Stear  1976). 
In  EKF  framework,  6  can  be  treated  as  a  random  process  with  the  following  model: 

0(k+i)  =  e(k)  +  cM,  (13) 


where  sd(k)  is  a  zero-mean  Gaussian  process  with  covariance  matrix  Sd. 

Table  1  summarizes  the  joint  estimation  of  pharmacokinetic  rates  and  ICG  concentration 
in  different  compartments.  In  table  1,  Cd(k\k  —  1)  is  the  state  estimate  propagation  at  step 
k  given  all  the  measurements  up  to  step  k  —  1,  C d(k)  is  the  state  estimate  update  at  step 
k,  Pjt.i-i  denotes  the  error  covariance  propagation  at  step  k  given  all  the  measurements  up  to 
step  k  —  1 ,  is  the  error  covariance  update  at  step  k,  Sd  is  the  preassigned  covariance  matrix 
of  sd(k),  J/t  is  the  Jacobian  matrix  due  to  iterative  linearization  of  the  state  equation  at  step 
k,  G/c  is  the  recursive  Kalman  gain  at  step  k,  R,/  is  the  covariance  matrix  of  the  measurements, 
Qd  is  the  covariance  matrix  of  the  concentration  vector  and  I  is  the  identity  matrix.  A  detailed 
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Figure  6.  Pharmacokinetic  rate  images,  (a)  km  and  (b)  k0 ut  for  case  1.  The  k[n  images  are  shown 
with  approximate  tumor  location  and  size. 


Table  1.  EKF  algorithm  for  simultaneous  estimation  of  states  and  parameters. 
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discussion  of  the  extended  Kalman  filtering  algorithm,  and  the  initialization  of  the  parameters, 
concentrations  and  covariance  matrices  can  be  found  in  Alacam  et  al  (2006). 


4.  Spatially  resolved  ICG  pharmacokinetic  rate  analysis  of  in  vivo  breast  data 

4.1.  Apparatus 

In  this  work,  we  use  the  data  collected  with  a  continuous  wave  (CW)  NIR  imaging  apparatus. 
The  apparatus  has  16  light  sources,  which  are  tungsten  bulbs  with  less  than  1  W  of  output 
power.  They  are  located  on  a  circular  holder  at  an  equal  distance  from  each  other  with  22.5° 
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Figure  7.  Pharmacokinetic-rate  images,  (a)  km  and  (b)  kou t  for  case  2.  The  km  images  are  shown 
with  approximate  tumor  location  and  size. 


X  [cm]  x  M 

(a)  (b) 


Figure  8.  Pharmacokinetic-rate  images,  (a)  km  and  (b)  k0 ut  for  case  3.  The  A:in  images  are  shown 
with  approximate  tumor  location  and  size. 


apart.  Sixteen  detectors,  namely  silicon  photodiodes,  are  situated  in  the  same  plane.  The 
breast  is  arranged  in  a  pendular  geometry  with  the  source-detector  probes  gently  touching  its 
surface.  Figure  5  illustrates  the  configuration  of  the  apparatus  and  the  configuration  of  the 
detectors  and  the  sources  in  a  circular  plane.  Note  that  sources  and  detectors  are  co-located. 
The  detectors  use  the  same  positions  as  the  sources  to  collect  the  light  originating  from  one 
source  at  a  time.  Only  the  signals  from  the  farthest  11  detectors  are  used  in  the  analysis. 
For  example,  when  source  1  is  on,  the  data  are  collected  using  detectors  4-14.  This  provides 
sufficient  number  of  source-detector  readings  (176  readings)  to  reconstruct  Afia  images  at 
each  time  instant.  A  band  pass  filter  at  805  nm,  the  absorption  peak  of  ICG,  is  placed  in  front 
of  the  sources  to  select  the  desired  wavelength.  A  set  of  data  for  one  source  is  collected  every 
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Figure  9.  ICG  concentration  images  in  plasma  for  case  1  for  (a)  246.4th,  (b)  334.4th  and 
(c)  422.4th  s. 


~500  ms.  The  total  time  for  a  whole  scan  of  the  breast  including  16  sources  and  16  detectors 
is  ~8.8  s.  A  more  detailed  explanation  of  the  apparatus  and  the  data  collection  procedure  can 
be  found  in  (Nioka  el  al  1997). 

4.2.  Tumor  information  and  protocol 

Three  different  patients  with  different  tumor  types  are  included  in  this  study.  Measurements 
are  made  before  the  biopsy  to  avoid  modification  of  the  blood  volume  and  flow  in  the  tumor 
region.  First  case  (case  1)  is  fibroadenoma,  which  corresponds  to  a  mass  estimated  to  be 
1-2  cm  in  diameter  within  a  breast  of  9  cm  diameter  located  at  6-7  o’clock.  Second  case 
(case  2)  is  adenocarcinoma  corresponding  to  a  tumor  estimated  to  be  2-3  cm  in  diameter 
within  a  breast  of  7.7  cm  diameter  located  at  4-5  o’clock.  The  third  case  (case  3)  is  invasive 
ductal  carcinoma,  which  corresponds  to  a  mass  estimated  to  be  4  cm  by  3  cm  located  at 
6  o’clock.  Table  2  describes  the  tumor  information  for  each  patient.  A  priori  information  on 
the  location  and  size  of  the  tumor  was  obtained  by  palpation  and  the  diagnostic  information 


848 


B  Alacam  et  al 


4  2  0  -2  -4 

X  [cm] 


(c) 

Figure  10.  ICG  concentration  images  in  the  EES  for  case  1  for  (a)  246.4th,  (b)  334.4th  and 
(c)  422.4th  s. 


Table  2.  Tumor  information  for  each  patient. 


Tumor  type 

Tumor  size 

Tumor  location 

Case  1 

Fibroadenoma 

1-2  cm 

6-7  o ‘clock 

Case  2 

Adenocarcinoma 

2-3  cm 

4-5  o ‘clock 

Case  3 

Invasive  ductal  carcinoma 

4  cm  by  3  cm 

6  o ‘clock 

was  derived  a  posteriori  from  biopsy  and  surgery.  ICG  is  injected  intravenously  by  bolus  with 
a  concentration  of  0.25  mg  kg1  of  body  weight.  Data  acquisition  started  before  the  injection 
of  ICG  and  continued  for  10  min. 

4.3.  Results  and  discussion 

Using  the  CW  imager  described  above,  source-detector  readings  were  collected  from  different 
angles  for  each  patient.  Differential  absorption  coefficient  images  were  reconstructed  based 
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Figure  11.  ICG  concentration  images  in  plasma  for  case  2  for  (a)  228.8th,  (b)  316.8th  and 
(c)  404.8th  s. 


on  DDOT  forward  model  given  in  equations  (1 )— (4)  with  to  set  to  zero.  Using  the  linear 
relationship  (5)  between  ICG  concentration  and  absorption  coefficient,  ICG  concentration 
images  were  obtained  for  each  case.  A  sample  set  of  ICG  concentration  images  for  the 
selected  time  instants  are  shown  in  figures  1-3  for  cases  1-3,  respectively.  Although  only  nine 
images  are  displayed,  there  are  approximately  50  images  for  each  case,  each  corresponding 
to  a  different  time  instant.  Each  image  is  composed  of  649  voxels.  Note  that  the  ICG 
concentration  images  in  figures  1-3  represent  the  bulk  ICG  concentrations  in  the  tissue,  not 
the  ICG  concentrations  in  plasma  or  the  EES  compartments. 

We  next  extracted  the  time  course  of  ICG  concentration  for  each  voxel.  As  an  example, 
figure  4  shows  the  time  course  of  ICG  concentrations  for  all  three  cases  for  a  specific  voxel 
in  the  tumor  region  (65th,  276th,  188th  voxel  for  cases  1,  2  and  3,  respectively).  We  then 
fit  the  two-compartment  model  to  each  time  course  data  using  the  EKF  framework;  and 
estimated  km,  kout,  ke im,  and  the  ICG  concentrations  in  plasma  and  the  EES.  We  chose  initial 
values  within  the  biological  limits  that  lead  to  minimum  norm  error  covariance  matrix.  The 
images  of  km  and  kout  for  each  case  are  shown  in  figures  6(a)-(b),  and  7(a)-(b),  8(a)-(b), 
respectively.  Additionally,  we  constructed  the  ICG  concentration  images  for  plasma  and 


850 


B  Alacam  et  al 


(.iM 


-3-2-1  0  1  2  3 

X  [cm] 


(a) 


(iM 


-3-2-1  0  1  2  3 


X  [cm] 


(b) 


Figure  12.  ICG  concentration  images  in  the  EES  for  case  2  for  (a)  228.8th,  (b)  316.8th  and 
(c)  404.8th  s. 


the  EES  compartments.  Figures  9-14  show  the  ICG  concentration  in  plasma  and  the  EES 
for  three  different  time  instants  for  cases  1,  2  and  3,  respectively.  Our  results  show  that 
the  pharmacokinetic  rates  are  higher  around  the  tumor  region  agreeing  with  the  fact  that 
permeability  increases  around  the  tumor  region  due  to  compromised  capillaries  of  tumor 
vessels.  We  also  observed  that  ICG  concentrations  in  plasma  and  the  EES  compartments  are 
higher  around  the  tumors  agreeing  with  the  hypothesis  that  around  the  tumor  region  ICG  may 
act  as  a  diffusible  extravascular  flow  in  leaky  capillary  of  tumor  vessels. 

Using  the  a  priori  and  a  posteriori  information  on  the  location,  and  the  size  of  the  tumors, 
we  plotted  an  ellipse  (or  a  circle)  to  identify  the  approximate  location  and  size  of  the  tumor 
in  the  pharmacokinetic-rate  images.  We  note  that  the  radii  of  the  ellipses  were  chosen  large 
enough  to  include  the  tumor  boundaries.  Figures  6(a),  7(a)  and  8(a)  present  the  km  images 
with  approximate  tumor  location  and  size  for  cases  1,  2  and  3,  respectively.  The  consistency 
of  the  bright  regions  in  the  km  images,  and  circular/elliptical  regions  drawn  based  on  the  a 
priori  and  a  posteriori  information  shows  that  the  pharmacokinetic-rate  images  may  provide 
good  localization  of  tumors. 
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Figure  13.  ICG  concentration  images  in  plasma  for  case  3  for  (a)  246.4th,  (b)  378.4th  and 
(c)  5 10.4th  s. 


Table  3.  Mean  and  standard  deviation  of  pharmacokinetic  rates  for  the  tumor  region  and  outside 
the  tumor  region. 


km  (sec 

1 10  2) 

^out  (sec 

1 10~2) 

^elm  (sec 

‘lO^3) 

Inside 

Outside 

Inside 

Outside 

Inside 

Outside 

Case  1 

2.14  ±0.018 

0.73  ±0.011 

1.24  ±0.069 

0.43  ±  0.013 

4.11  ±  0.057 

3.87  ±0.012 

Case  2 

2.92  ±  0.076 

1.14  ±0.052 

1.58  ±0.051 

0.65  ±  0.036 

3.94  ±0.081 

4.12  ±0.047 

Case  3 

6.87  ±  0.093 

3.06  ±0.015 

4.96  ±  0.048 

1.66  ±0.072 

4.49  ±  0.056 

4.46  ±  0.081 

The  histograms  of  km  and  kout  images  for  the  tumor  region  (as  indicated  by 
circular/elliptical  regions)  and  outside  the  tumor  region  are  shown  in  figures  15(a)-(c)  and 
figures  16(a)-(c),  respectively.  Note  that  all  nonzero  voxels  outside  the  elliptical  region 
constitute  ‘outside  the  tumor  region’.  The  solid  curves  in  figures  15  and  16  show  the  Gaussian 
fit.  The  histograms  and  their  Gaussian  fits  in  figures  15  and  16  show  that  the  mean  and  the 
standard  deviation  of  km  and  lcout  values  are  different  for  the  tumor  and  outside  the  tumor  region. 
Table  3  tabulates  the  mean  values  (±  spatial  standard  deviation)  of  the  pharmacokinetic  rates 
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Figure  14.  ICG  concentration  images  in  the  EES  for  case  3  for  (a)  246.4th,  (b)  378.4th  and 
(c)  5 10.4th  s. 


Table  4.  Bulk  pharmacokinetic  rates  extracted  from  the  entire  breast  tissue. 


kjn  (sec  1 10  2) 

*out  (sec  1 10  2) 

^elm  (sec  *10  3) 

Case  1 

0.84  ±0.013 

0.62  ±0.017 

3.66  ±  0.042 

Case  2 

2.01  ±  0.022 

0.83  ±0.012 

4.01  ±  0.054 

Case  3 

4.06  ±  0.072 

3.36  ±0.051 

4.37  ±  0.052 

for  the  tumor  region  and  outside  the  tumor  region  for  all  three  cases.  The  pharmacokinetic 
rates  are  higher  for  case  3  (invasive  ductal  carcinoma),  for  both  the  tumor  region  and  outside 
the  tumor  region  as  compared  to  case  2  (adenocarcinoma).  Similarly,  the  kinetic  rates  are 
higher  for  case  2  (adenocarcinoma),  as  compared  to  case  1  (fibroadenoma)  for  both  the  tumor 
region  and  outside  the  tumor  region.  This  observation  shows  that  high  mean  values  of  km  and 
koul  may  be  indicative  of  tumor  aggressiveness. 

To  understand  the  value  of  pharmacokinetic  rate  imaging  as  compared  to  the  bulk 
pharmacokinetic  rate  analysis,  we  averaged  the  concentration  images  spatially,  and  obtained 
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Figure  15.  The  histograms  of  k[n  for  (a)  case  1,  (b)  case  2,  (c)  case  3  for  the  tumor  region  (light 
gray)  and  outside  (blue/dark  gray)  the  tumor  region  (as  indicated  by  circular/elliptical  regions). 
The  solid  lines  in  figures  show  the  Gaussian  fit. 


a  bulk  concentration  value  for  each  time  instant.  We  then  formed  a  time  curve  for  the  bulk 
ICG  concentrations.  Next,  we  fit  the  two-compartment  model  to  the  resulting  time  curves  and 
estimated  the  bulk  pharmacokinetic  rates.  Table  4  tabulates  the  bulk  pharmacokinetic  rates 
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Figure  16.  The  histograms  of  £0ut  for  (a)  case  1,  (b)  case  2,  (c)  case  3  for  the  tumor  region  (light 
gray)  and  outside  (blue/dark  gray)  the  tumor  region  (as  indicated  by  circular/elliptical  regions). 
The  solid  lines  in  figures  show  the  Gaussian  fit. 


for  each  patient.  To  compare  the  bulk  rates  with  spatially  resolved  rates,  in  figures  17  and 
18,  the  bulk  pharmacokinetic  rates  are  overlaid  on  the  histograms  of  the  pharmacokinetic  rate 
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Figure  17.  Solid  lines  (blue)  shows  bulk  kin  rates  for  (a)  case  1,  (b)  case  2,  (c)  case  3  together  with 
the  histogram  fits.  The  dashed  (red)  line  indicates  the  Bayesian  minimum  error  classifier  threshold. 

images.  The  dotted  line  shows  the  Bayesian  minimum  error  classifier  threshold  (the  value 
corresponding  to  the  intersection  of  the  histograms)  (Fukunaga  1990)  for  each  case.  We  see 
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Figure  18.  Solid  lines  (blue)  shows  bulk  kout  rates  for  (a)  case  1 ,  (b)  case  2,  (c)  case  3  together 
with  the  histogram  fits.  The  dashed  (red)  line  indicates  the  Bayesian  minimum  error  classifier 
threshold. 


that  for  case  1,  the  bulk  rates  of  and  koui  are  both  classified  as  healthy  tissue  (outside  the 
tumor  region).  For  case  2,  km  is  classified  as  cancerous  tissue  (in  the  tumor  region)  and  kout  is 
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classified  as  healthy  tissue.  Similarly  for  case  3,  km  is  classified  as  healthy  and  k0llt  is  classified 
as  cancerous  tissue.  This  indicates  that  spatially  resolved  rates  may  provide  more  consistent 
and  superior  information  than  the  bulk  rates. 

5.  Conclusion 

In  this  study,  we  presented  a  method  of  forming  pharmacokinetic-rate  images  and  reported 
pharmacokinetic  rate  images  of  ICG  for  three  patients  with  breast  tumors.  To  form 
pharmacokinetic  rate  images,  we  first  obtained  a  sequence  of  ICG  concentration  images  using 
the  differential  diffuse  optical  tomography  technique.  We  next  employed  the  two-compartment 
model,  and  estimated  the  pharmacokinetic  rates  and  concentrations  in  each  compartment  for 
each  voxel  using  the  EKF  framework.  We  have  shown  in  our  prior  work  (Alacam  et  al  2006) 
that  the  EKF  framework  has  a  number  of  advantages  in  pharmacokinetic  rate  estimation,  some 
of  which  include  robust  estimation  in  the  presence  of  measurement  noise  and  dynamic  model 
uncertainties. 

We  formed  the  pharmacokinetic  rate  images  using  the  in  vivo  data  obtained  from  three 
patients  with  breast  tumors.  We  also  obtained  bulk  pharmacokinetic  rates  for  each  patient. 
Both  spatially  resolved  and  bulk  rates  show  that  high  values  of  km  and  kout  may  be  indicative 
of  tumor  aggressiveness.  Along  with  the  pharmacokinetic  rates,  we  also  estimated  the  ICG 
concentrations  in  plasma  and  EES  compartments.  We  observed  that  ICG  concentrations  in 
plasma  and  the  EES  compartments  are  higher  in  the  tumor  region  agreeing  with  the  hypothesis 
that  around  the  tumor  region  ICG  may  act  as  a  diffusible  extravascular  flow  in  leaky  capillary 
of  tumor  vessels. 

Comparison  of  spatially  resolved  and  bulk  ICG  pharmacokinetic  rates  show  that  ICG 
pharmacokinetic  imaging  may  provide  more  consistent  and  superior  information  than  bulk 
ICG  pharmacokinetic  rates. 

While  the  available  patient  data  are  limited  to  perform  a  full  scale  receiver  operating 
characteristic  study,  clearly,  pharmacokinetic  rate  imaging  provides  a  new  tool  to  investigate 
and  improve  breast  cancer  diagnosis,  staging,  and  treatment  monitoring.  This  includes 
extraction  of  new  quantitative  features  from  ICG  pharmacokinetic  rate  images,  within  patient 
comparison  of  these  features,  and  statistical  analysis  of  spatial  distribution  of  pharmacokinetic 
rates.  We  leave  for  future  work  to  collect  sufficient  number  of  patient  data,  and  to  fully 
investigate  the  value  of  ICG  pharmacokinetic  rate  imaging  for  breast  cancer  diagnosis,  staging, 
and  treatment  monitoring. 
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